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1.  INTRODUCTION 


In  the  first  year  of  this  study  we  showed  that  the  dyadic  transform  was  not  sufficient  for  our  case  mass 
of  detection  [1].  Last  year,  we  developed  multi-scale  adaptive  histogram  equalization  (MSAHE)  that 
achieved  a  global  contrast  enhancement  by  adjusting  contrast  locally  through  reassigning  a  central 
pixel  the  value  through  a  local  histogram  equalization  mapping  function.  This  past  six  months  we 
have  implemented  a  level  set  method  of  segmentation  for  the  detection  of  masses  within  a  multi-scale 
expansion.  This  report  provides  an  overview  of  this  method  and  describes  progress  to  date.  Since  this 
is  a  final  report,  we  will  summarize  previous  accomplishments  and  status  of  the  project  overall,  at  the 
end. 

Traditional  methods  of  segmentation,  such  as  pixel-based  clustering,  region  growing,  and  edge 
detection,  requires  additional  pre-processing  and  post-processing  as  well  as  a  considerable  amounts  of 
expert  intervention  or  information  of  the  objects  of  interest.  Furthermore  the  subsequent  analysis  of 
segmented  objects  is  hampered  by  primitive,  pixel  or  voxel  level  representations  from  region-based 
segmentation  [1]. 

Deformable  models,  on  the  other  hand,  provide  an  explicit  representation  of  the  boundary  and  the 
shape  of  an  object.  They  combine  several  desirable  features  such  as  inherent  connectivity  and 
smoothness,  which  counteract  noise  and  boundary  irregularities,  as  well  as  tbe  ability  to  incorporate 
knowledge  about  the  object  of  interest  [1,3]  [4].  However,  parametric  deformable  models  have  two 
main  limitations.  First,  in  situations  where  the  initial  model  and  desired  object  boundary  differ  greatly 
in  size  and  shape,  the  model  must  be  re-parameterized  dynamically  to  faithfully  recover  the  object 
boundary.  The  second  hmitation  is  that  it  has  difficulty  dealing  with  topological  adaptation  such  as 
splitting  or  merging  model  parts,  a  useful  property  for  recovering  either  multiple  objects  or  an  object 
with  unknown  topology.  This  difficulty  is  caused  by  the  fact  that  a  new  parameterization  must  be 
constructed  whenever  the  topology  change  occurs,  which  requires  sophisticated  schemes  [5,  6]. 

In  the  body  of  this  report  an  alternative  approach  of  coping  with  weaknesses  in  existing  methods  of 
segmentation  of  masses  is  described. 


2.  BODY 


DETECION  OF  MASSES  VIA  EVOLUTION  OF  LEVEL  SET  BOUDARIES 


In  existing  level-set  methods,  the  gradient  information  is  used  as  a  stopping  criterion  for  curve 
evolution,  and  also  provides  the  attracting  force  to  the  zero  level-set  from  the  target  boundary. 
However,  in  a  discrete  implementation,  the  gradient-based  term  can  never  fully  stop  the  level-set 
evolution  even  for  ideal  edges,  leakage  is  often  unavoidable.  Also  the  effective  distance  of  the 
attracting  force  and  blurring  of  edges  become  a  trade-off  in  choosing  the  shape  and  support  of  the 
smoothing  filter.  The  proposed  homogeneity  measurement  provides  easier  and  more  robust  edge 
estimation,  and  the  possibility  of  fully  stopping  the  level-set  evolution.  The  homogeneity  term 
decreases  from  a  homogenous  region  to  the  boundary,  which  dramatically  increases  the  effective 
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distance  of  the  attracting  force  and  also  provides  an  additional  measurement  of  the  overall 
approximation  to  the  target  mass  boundary.  Therefore,  it  provides  a  reliable  criterion  of  adaptively 
changing  the  advent  speed.  By  using  this  term,  the  leakage  problem  was  avoided  effectively  in  most 
cases  compared  to  traditional  level-set  methods.  The  computation  of  the  homogeneity  operator  is  fast 
and  can  be  done  within  seconds  on  a  PC  workstation. 


Curvature  evolution  and  Level  set 

Level  set  segmentation  [7,  8],  also  referred  as  geometric  deformable  models,  provides  an  elegant 
solution  to  address  the  primary  limitations  of  parametric  deformable  models.  In  the  2D  case,  the 
boimdary  of  an  object  is  implicitly  represented  as  the  zero-level  set  of  a  time  dependent  2D  function, 
which  is  usually  called  the  level  set  function.  A  useful  property  of  this  approach  is  that  the  level  set 
function  remains  a  valid  function  while  the  embedded  curve  (the  zero  level  set)  can  change  its 
topology. 

The  evolution  equation  for  the  level  set  function  ^(A:,y,0  takes  on  the  following  formula  [9]: 

^+F|V«(|  =  0  (1) 

The  evolution  of  the  level  set  function  was  determined  by  the  speed  function  F.  As  an  example, 
imagine  that  given  an  initial  closed  curve  that  is  evolving  under  three  simultaneous  motions.  First,  it  is 
expanding  with  a  constant  speed  in  its  normal  direction;  second,  it  is  moving  witii  a  speed  proportional 
to  its  curvature;  third,  it  is  being  passively  advected  by  an  imderlying  velocity  field  whose  direction 
and  strength  depend  on  position  and  time,  but  not  on  the  front  itself  This  entire  motion  can  then  be 
written  in  terms  of  the  speed  function  as  an  explicit  level  set  scheme; 

P  =  Pprop  +  ^curv  +  ^adv  (2) 

where  =  F^  is  the  propagation  expansion  speed,  F^^  =  -fir  is  the  dependence  of  the  speed  on  the 

curvature  k,  and  F^^^  =  U{x,y,tynis  the  advection  speed,  where  nis  the  normal  to  the  front.  The 
PDE  in  (1)  can  be  solved  with  entropy-satisfying  schemes  given  the  speed  function. 

A  small  modification  version  of  (2)  gives  the  general  formula  for  level  set  segmentation: 

(l>,+gj{\-SK)\V(l>\-INP-V<l>  =  0  (3) 


In  the  term  \-£k  ,  the  uniform  expansion  with  speed  1  corresponds  to  the  inflation  force  used  by 
Cohen  [10].  The  diffusive  term  f/c  smoothes  out  the  high  curvature  regions  and  has  the  same 
regularizing  effect  as  the  internal  deformation  energy  term  in  thin-plate-membrane  splines  [2].  The 
term  g^was  computed  from  the  image  data,  and  provide  a  halting  criterion  for  the  speed  function,  the 
value  of  gj  should  be  between  0  and  1,  and  ideally,  with  0  on  the  boundary  and  1  within  the 
homogeneous  region  (either  within  or  outside  the  object).  Typically,  it  can  be  estimated  by  the 
gradient; 


giix,y)  = 


_ 1 _ 

l  +  |V(G,*7(x,y)) 


(4) 


where  the  expression  G^*/ denotes  the  image  convolved  with  a  Gaussian  smoothing  filter  whose 
characteristic  width  is  cr .  The  term  V(G^  *I{x,y))  is  essentially  zero  except  where  the  image  gradient 
changes  rapidly,  in  which  case  the  value  becomes  large.  Thus,  gj  is  close  to  unity  away  from 
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boimdaries  and  drops  to  zero  near  sharp  changes  in  the  image  gradient.  These  changes  presumably 
correspond  to  the  edges  of  the  desired  shape.  In  other  words,  the  first  term  in  (3)  anticipates  steep 
drops  in  the  image  gradient,  and  retards  the  evolving  fi-ont  fi*om  passing  out  of  file  desired  region.  The 
second  term  in  (3)  is  a  force  which  attracts  the  surface  towards  the  boundary,  which  has  a  stabilizing 
effect,  especially  when  there  is  a  large  variation  in  the  image  gradient  value.  This  term  denotes  the 
projection  of  an  attractive  force  vector  on  the  surface  normal.  This  force,  introduced  in  [1 1],  is  realized 
as  the  gradient  of  a  potential  field.  Here; 

/>  =  -|V(G,./)|  (5) 

attracts  the  surface  to  the  edges  in  the  image,  the  coefficient  controls  the  strength  of  this  attraction. 

Apparently,  both  terms  depend  on  the  edge  map  V(G^  *1) ,  and  the  quality  of  this  edge  estimation 
determines  performance  of  the  segmentation. 

‘^Scale  Map”  based  on  homogeneity  measurement 

“Scale”  is  a  fimdamental,  well-established  concept  in  image  processing  [12,  13].  The  premise  behind 
this  concept  is  to  consider  the  local  size  of  the  object  in  carrying  out  whatever  local  operations  that  is 
to  be  carried  out  on  the  image.  It  has  previously  been  used  as  a  metric  of  local  homogeneity  [13]. 

“Object  Scale”  in  an  image  C  at  any  pixel  c  was  defined  as  the  radius  r(c)  of  the  largest  hyper  ball 
centered  at  c  which  lies  entirely  in  the  object  region  [14]. 


A  hyper  ball  (c)  centered  at  pixel  c  is  a  collection  of  pixels  around  c,  i.e.  (c)  =  e  C  |  ||c  -  e||  <  r| . 

For  a  ball  of  any  radius  k  centered  at  c,  a  fi-action  fimction  FO^{c)  was  defined,  which  indicates 

whether  the  fraction  of  the  ball  boundary  is  sufficiently  homogeneous  to  the  inside  region  of  the  ball: 


FO,(c)  = 


where  \B^{c)-B^_^(c)\is  the  number  of  pixels  in  B^(c)-B^_^(c)and  W  is  the  homogeneity  fimction, 

which  measures  the  similarity  of  two  pixels  based  on  their  pixel  value  in  the  image  f(x).  Some 
typically  used  homogeneity  fimctions  for  this  sake  are  illustrated  in  Figure  1. 


Figure  1;  Examples  of  homogeneity  functions  for  computing  the  fraction  fimction  in  (6). 
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The  mathematical  formulas  for  the  homogeneity  functions  in  Figure  1  are: 


(a)  W(x)  = 


Jl,  0<x<a 
[0,  x>a 


(b)  W(x)  = 


1, 

b-x 


b-a 

0, 


0<  a:<  a 
a<x<b 
x>b 


(c)  W(x)  =  e^^^\k>0 


(7) 


By  setting  a  fixed  thresholdO<t  <1,  the  scale  map  at  a  pixel  c  can  be  computed  by  the  following 
pseudo-code: 

begin 

k=l; 

while  FOj^{c)>t  do 

k=k+l; 

endwhile 

r(c)=k; 


Methodology  for  computing  an  edge  map  based  on  homogeneity  metrics 

The  scale  map  described  above  provides  a  robust  homogeneity  measurement  by  incorporating  a 
tolerance  level  t.  For  example,  if  we  use  r  =  J^,ina3x3  neighborhood  of  a  pixel  c  in  a  2D  scene,  we 

allow  one  out  of  the  eight  neighboring  pixels  to  belong  to  a  different  object  (to  account  for  noise)  but 
still  consider  the  neighborhood  to  be  entirely  within  the  same  object.  This  actually  provides  a 
mechanism  of  denoising  within  the  homogeneity  measurement. 

An  edge  was  defined  as  a  region  of  an  image  in  which  the  pixel  value  changes  significantly  over  a 
short  distance  [15].  Therefore  the  edge  represents  a  region  that  has  significant  lower  homogeneity. 
Thus  a  homogeneity  measurement  certainly  provides  edge  information  within  an  image. 

Figure  2  shows  the  “scale  map”  computed  by  the  algorithm  above  using  different  parameters. 
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(a)  (b) 

Figure  2:  scale  map  of  breast  radiograph  image,  (a)  original  mammogram  (256  gray  levels), 
(b)  scale  map  computed  with  parameter,  =500,  t=0.8. 


Here  we  use  equation  7(c)  as  the  homogeneity  function  for  computation.  The  two  parameters  for 
computing  the  scale  map  are  the  shape  parameter  of  the  homogeneity  function  (A),  and  die  tolerance 
value  for  homogeneity  measurement  (r).  Parameter  k  determines  how  much  variation  in  the  pixel  value 
is  tolerated  in  terms  of  homogeneity.  Paramter  t,  as  discussed  previously,  determines  how  much  noise 
we  want  to  ignore.  Figure  2(b)  shows  the  appearance  of  an  edge  map  with  lower  values  on  the 
boundary  and  higher  values  on  the  homogeneous  regions.  As  shown  in  Figure  3,  when  we  scale  the 
pixel  values  of  this  edge  map  to  the  interval  [0,1],  it  can  be  used  effectively  as  an  image-based  term  for 
level  set  segmentation  of  masses. 

Level  set  segmentation  using  a  homogeneity  edge  map 

The  level  set  evolution  function  (3),  without  the  attracting  speed  term,  can  be  written  as: 

?>,+«, =  (8) 
where  the  first  term  in  (8)  provide  the  expansion  speed  along  the  normal  direction  of  the  curve,  the 
magnitude  is  gj ,  therefore,  in  the  non-boundary  region,  we  want  the  edge  map  gj  to  be  as  large  as 
possible,  so  that  the  curve  evolution  can  converge  to  the  real  boundary  quickly.  Ideally,  gj  will  be  zero 
on  the  estimated  boundary,  so  that  when  the  ciuwe  reaches  it,  the  evolution  function  becomes^,  =  0 , 
and  reaches  equilibrium.  But  since  the  gradient  V(G^*/)  will  never  be  infinity  in  the  discrete 
implementation,  gj  will  never  become  exactly  zero.  Thus  care  and  reliable  methods  are  needed  to  stop 
the  level  set  evolution  when  it  reaches  the  estimated  boundary. 

In  this  research,  we  used  ttie  edge  map  defined  by  a  homogeneity  metric,  where  the  value  of  the  edge 
map  can  be  zero  on  a  sharp  boimdary.  This  provided  a  more  reliable  stopping  criterion  than  traditional 
gradient  operators  alone.  However,  because  of  noise  tolerance,  a  weak  boundary  could  be  non-zero. 
Also,  missing  boundaries  may  occur  due  to  the  angle  of  an  X-ray  projection,  therefore,  boimdary 
leakage  may  happen  if  the  boundary  definition  in  the  original  image  is  not  perfectly  clear.  This  in  part 
motivated  our  development  of  the  MSAHE  algorithm  described  in  the  beginning  sections  in  this  report. 
The  attraction  term  (5)  introduced  in  [1 1]  is  used  to  pull  back  the  curve  when  it  passes  a  boundary.  We 
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added  an  adaptive  indication  term  that  shuts  down  the  expansion  speed  when  the  curve  became  close  to 
the  estimated  boundary.  The  evolution  function  we  used  was: 

<l>,+gj-5-\V<l>\-SK\V<l\-INP*V(l>  =  0.  (9) 

The  reason  why  we  took  out  gj  from  the  second  term  is  that  when  the  curve  approaches  the  boundary, 
the  value  of  became  very  small,  and  therefore  the  smoothing  effect  was  eliminated.  When  using  the 
attracting  term,  the  cmve  always  appeared  noisy  due  to  imperfect  boundaries  in  the  edge  map.  This 
strategy  is  analogous  to  those  used  in  parametric  deformable  models  that  always  keep  a  constant 
weighted  elastic  internal  force. 

5  is  an  global  indication  function  such  that  when  the  zero  level  set  of  is  close  enough  to  the 
estimated  boundary,  ^  =  0  and  otherwise  S  =  l.  Looking  at  the  computation  of  the  homogeneity  map, 
intuitively  the  pixel  value  in  the  edge  map  is  decreasing  from  the  homogeneous  region  to  the  boundary, 
and  reaches  the  minimum  at  the  boimdary.  Therefore  the  pixel  value  of  the  edge  map  actually  gives  the 
information  of  how  close  a  certain  pixel  is  to  the  boimdary.  By  averaging  the  values  in  the  edge  map 
over  the  pixels  that  represent  the  zero  level  set,  or  finding  the  maximum  value,  we  can  ascertain  how 
close  the  whole  curve  is  to  the  target  boundary.  A  threshold  is  chosen  so  that  when  this  value  is  smaller 
than  some  threshold  =  0 ,  and  otherwise  <5  =  1. 
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3.  KEY  RESEARH  ACCOMPLISHMENTS 


During  the  past  six  months  of  the  project  we  incorporated  an  expansion  based  method  of  detecting 
masses  in  digital  radiographs. 

•  We  implemented  a  level-set  method  of  segmentation  that  made  use  of  a  local  homogeneity 
operator  for  the  detection  of  subtle  masses  in  digital  mammography. 

•  This  expansion  based  method  is  integrated  into  our  previously  described  multi-scale  expansion 
framework  and  will  be  tested  using  a  local  database  of  digital  mammograms  as  part  of  a 
planned  validation  study. 


In  our  proposal  we  made  the  conjecture  that  existing  methods  of  analysis  computed  (efficiently)  on 
dyadic  scales  were  not  sufficient  for  the  detection  of  masses:  a  lesion  may  be  too  blurry  at  one  scale, 
and  too  precise  at  the  next  finer  (dyadic)  scale.  In  the  first  part  of  this  study,  we  answer  the  question  “Is 
it  sufficient  to  work  with  dyadic  scales,  or  is  there  an  absolute  need  to  compute  coefficients  between 
the  scales?”  Indeed,  during  the  course  of  this  investigation  we  have  clearly  shown  that  there  is  a 
significant  advantage  in  the  capturing  the  morphology  of  arbitrary  sized  masses  when  using  finer 
“grained”  expansions. 

Throughout  this  study  we  have  avoided  development  of  not  “reinventing  the  wheel”  whenever 
possible.  We  choose  to  use  or  modify  existing  libraries  and  programs  readily  available  within  the 
mathematical  community.  For  example,  we  modified  under  Matlab  several  existing  LastWave 
algorithms,  including  the  Discrete  Wavelet  Transform  in  two  dimensions  without  downsampling,  using 
the  Algorithme  ^  Trous  algorithm.  An  ancillary  benefit  to  this  approach  is  that  when  this  code  is  made 
available  to  the  research  community  it  will  be  easy  to  use  having  been  built  upon  “freeware”  and 
commercially  available  programming  environments.  All  of  the  programs  for  computing  the  expansion 
and  detection  algorithms  (written  in  “C”  and/or  MatLab  code  )  during  this  the  course  of  the  study  are 
available  upon  request  through  our  web  site:  bil.bme.colimibia.edu”. 

We  compared  in  one  dimension  the  CWT  and  the  DWT  in  order  to  show  a  proof  of  concept  concerning 
any  advantage  of  pursuing  refinement  of  scale.  We  processed  phantom  masses,  and  ID  intensity 
profiles  of  real  masses  mammograms  to  evaluate  feasibility.  In  order  to  identify  the  best  scale,  we 
evaluated  the  use  of  maxima  of  die  coefficients  and  a  correlated  model  using  three  masses  of  different 
size.  We  then  evaluated  the  shape  of  the  “Mexican  hat”  for  suitability  in  a  matched  filtering  detection 
paradigm.  During  the  third  year  we  applied  enhancement  algorithms  as  a  preprocessing  step  and 
introduced  the  notion  of  “voices”  which  allowed  us  to  compute  representations  of  masses  in  between 
octaves  of  the  expansion. 


4.  REPORTABLE  OUTCOMES 


The  following  publications  are  presented  as  reportable  outcomes  of  the  investigation.  For 
convenience,  we  have  included  selected  copies  of  these  manuscripts  along  with  this  report. 
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Publications: 

(0)  Y,  Jin,  A,  Laine  and  C.  Imielinska,  "An  adaptive  speed  term  based  on  homogeneity  for  level-set 
segmentation,"  Medical  Imaging,  Proceedings  of  SPIE.,  Vol.  4684  (1),  pp.  383-390,  Feb.  2002,  San 
Diego  CA. 

(1)  M.  A.  Birgen,  S.  Smith,  A.  F.  Laine,  “Detection  of  Masses  in  Mammography  Through  Redundant 
Exapansions  of  Scale,  ’’Proceedings  of  the  23rd  Annual  International  Conference  of  the  IEEE 
Engineering  in  Medicine  and  Biology  Society,  Istanbul,  Turkey,  October;  2001. 

(2)  R.  Mekle,  A.  Laine,  S.  Smith "  Evaluation  of  a  Multi-Scale  Enhancement  Protocol  for  Digital 
Mammography,"  Image-Processing  Techniques  For  Tumor  Detection,  R.  N.  Strickland,  Ed.,  Marcel 
Dekker,  New  York,  NY,  (2001)  pp.  155-186. 

(3) ,W.  Huda,  Y.  Jin  and  A.  Laine,  "Evaluation  of  Contrast  Enhancement  by  Digital  Equalization  in 
Digital  Mammography,"  World  Congress  on  Medical  Physics  and  Biomedical  Engineering,  Chicago, 
2000. 

(4)  R.  Mekle,  A.  Laine,  S.  Smith,  C.  Singer,  T.  Koenigsberg  and  M.  Brown,  "  Evaluation  of  a 
Multiscale  Enhancement  Protocol  for  Digital  Mammography,"  in  Wavelet  Applications  in  Signal  and 
Image  Processing  VUI,  A.  Aldroubi,  A.  F.  Laine,  M.  A.  Unser,  Eds.,  Proceedings  of  the  SPIE  Vol. 
4119,  pp.  1038-1049,  San  Diego,  2000. 

(5)  Koren,  A.  Laine,  S.  Smith,  E.  Nickoloff  and  F.  Taylor,  "Visualization  of  Memmography  via 
Fusion  of  Enhanced  Features,"  in  M.  Doi  (Editor),  First  International  Workshop  on  Computer-Aided 
Diagnosis,  Elsevier  Science,  Amsterdam,  1999,  pp.  287-303. 


List  of  Personnel: 

Andrew  Laine,  D.Sc.  (Principal  Investigator) 

Suzanne  Smith,  M.D.  (Director  of  Breast  Imaging) 

C.  Singer,  M.D.  (Radiologist,  Mammography  Reader) 

T.  Koenigsberg,  M.D.  (Radiologist,  Mammography  Reader) 
M.  Brown,  M.D.  (Radiologist,  Mammography  Reader) 
Edward  Nickolloff,  Sc.D.  (Chief  Medical  Physicist) 

Ralf  Mekle,  Ph.D.  (Graduate  Research  Assistant) 

Yinpeng  Jin,  M.S.  (Graduate  Research  Assistant) 

Mehmet.  A.  Birgen.  (Graduate  Research  Assistant) 
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5.  CONCLUSIONS 


The  major  problem  of  level  set  segmentation  is  boundary  leakage  when  weak  boundaries  or  occluded 
parts  of  a  mass  boundary  are  evident.  In  this  research,  we  kept  the  expansion  term  initially  for 
“pushing”  the  model  towards  the  boundary,  and  therefore  had  no  need  for  a  restricted  initialization. 
When  the  expansion  term  was  shut  down,  we  observed  that  it  prevented  boundary  leakage. 

The  definition  of  our  homogeneity  measurement  did  not  include  any  prior  knowledge  of  the 
dimensions,  and  therefore  can  be  easily  extended  to  higher  dimensions  and  it  is  computationally 
efficient  because  no  filtering/convolution  is  involved.  In  addition,  to  collect  all  the  pixels  around  the 
zero  level  set  needed  to  compute  the  indication  function  S  was  a  byproduct  of  the  procedure  for 
constructing  the  extension  speed  when  iterating  the  level  set  function.  Therefore  this  test  can  be  done 
whenever  the  construction  of  die  extension  speed  is  performed  without  extra  computation. 

Since  this  research  has  presented  a  new  edge  map,  other  methods  based  on  speed  terms  driven  by  edge 
maps  [16,  20]  can  be  designed  to  achieve  more  reliable  level  set  evolution  and  segmentation.  In  the 
future,  extensions  to  higher  dimensional  datasets  and  building  other  desirable  speed  terms  based  on 
“hyper  edge  maps”  and  other  information  of  the  underlining  image  are  possible. 

Overall  Summary 

The  idea  of  this  “Idea  Award”  was  to  detect  subtle  masses  in  mammograms  by  tuning  the  central 
frequency  and  width  of  a  basis  function  that  generates  overcomplete  expansions.  By  modeling  the 
shape  of  a  mass  through  this  flexibility  we  hoped  to  detect  small  and  subtle  masses  in  dense  breasts 
and  improve  the  chances  of  early  detection  in  screening  mammography.  In  the  first  part  of  our 
investigation,  we  evaluated  existing  tools  to  compute  overcomplete  expansions  of  multiscale  signals. 
We  compared  in  one  dimension  the  CWT  and  the  DWT  for  a  proof  of  concept  concerning  any 
advantage  of  pursuing  refinement  of  scale.  We  processed  phantom  masses,  and  ID  intensity  profiles 
of  real  masses  mammograms  to  evaluate  feasibility.  In  order  to  identify  the  best  scale,  we  evaluated 
the  use  of  maxima  singularities  and  a  correlated  model  using  three  masses  of  different  size.  Our  study 
answered  the  question  of  weather  of  not  dyadic  scales  were  sufficient  to  detect  masses  in  a  dense 
mammograms.  We  clearly  showed  that  reasonable  approximations  of  mass  shapes  could  be  obtained 
through  overcomplete  expansions  that  computed  voices  between  the  traditional  dyadic  scales. 

Our  study  of  one  dimension  cases  answered  the  question  of  weather  of  not  dyadic  scales  were 
sufficient  to  detect  masses  in  a  dense  mammograms.  We  showed  that  reasonable  approximations  of 
mass  shapes  could  be  obtained  through  overcomplete  expansions  of  a  continuous  wavelet  transform 
that  computed  voices  between  the  traditional  dyadic  scales. 

We  observed  in  mathematical  phantoms  and  real  masses  that  a  correlation  method  (between  a  model  of 
a  mass  and  the  values  of  the  computed  coefficients)  gave  approximately  the  same  results  when 
compared  to  the  maxima  method  (maximum  of  the  coefficients  at  each  scale).  We  developed  a  simple 
scheme  to  detect  masses  using  these  representations.  This  method  based  on  geometric  properties  of 
segmented  masses  within  each  expansion  was  shown  to  be  remarkably  stable. 
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ABSTRACT 

We  show  that  dyadic  scales  may  not  be  sufficient  for  the 
detection  of  masses  in  mammograms:  a  lesion  may  be  too 
blurred  on  one  scale,  and  then  too  fragmented  at  the  next. 
In  this  paper,  we  report  on  the  preliminary  evidence  of 
our  study  using  a  continuous  wavelet  transform  in  two 
dimensions  with  arbitrary  positioning  of  a  wavelet’s 
center  frequency  channel  tuned  to  the  mass  detection 
problem.  Our  goal  is  to  detect  masses  in  dense 
mammograms  whose  diameter  is  smaller  than  1  cm.  The 
aim  is  to  be  able  to  find  the  scale  where  the  mass  is  best 
represented  in  terms  of  analysis. 

1.  INTRODUCTION 

An  initial  study  in  one  dimension  helped  us  observe  that 
dyadic  scales  are  often  not  sufficient  to  detect  a  mass  in  a 
dense  mammogram  [3].  Below  we  show  this  by  a 
continuous  wavelet  transform,  which  computes  the 
decomposition  on  voices  between  traditional  dyadic 
scales. 

2.  METHODOLOGY 
2.1.  Voices  and  octave 

It  is  possible  to  expand  a  signal  more  finely  and  compute 
scales  between  octaves  of  traditional  dyadic  expansions 
by  voices  [1].  A  voice  constitutes  a  subdivision  of  an 
octave.  If  we  consider  a  wavelet  mother  y/ ,  a  family  of 

m 

wavelets  ii/^^(x)=a^^yAcC'x-rib^)  where  ao  is  the 
dilatation  parameter,  bo  is  die  translation  parameter, 
(m,  n)  G  are  possible.  In  the  dyadic  case,  ao=2  and 

m 

t>o=l)  ^^„(jc)  =  2  Decomposing  N  voices 

per  octave  means  creating  N  functions  y/n,m 
computing  the  frame 

•  (^3^)  ^  Analyzing  with  N 


voices  means  finding  N  different  frequency  channels, 
which  correspond  to  the  N  frequency  localizations  of 

...5  -0^  [2],  all  translated  by  the  same  step  (Fig.  lb). 
Such  a  lattice  can  be  viewed  as  the  superposition  of  N 
different  lattices  of  the  type  shown  in  Fig.  la,  stretched 
by  fixed  amounts  in  frequency.  For  example  a  possible 

choice  for  is  ^(x)=2"^ 

If  ,  which  we  assume  to  be  even,  peaks 

around  ±€0^^  then  |^?'“'  jwill  be  concentrated  around 

+2  ^  in  the  same  way  as  in  the  dyadic  case.  If 

/V 

'll)  has  two  peaks  in  frequency  at 
peaks  at  ±2^”  which  are  two  localization  centers  of 

^m,n- 

The  equation  computing  the  scale  for  source 
given  “octave”,  “current  voice”  and  “number  of  voices” 

octave+  current  _,voice 

IS  scale  =  2  number^voices  1-3]  ^ 

Moreover,  we  adopt  the  following  convention: 
the  first  octave  (octave  number  zero)  corresponds  to  the 

1 

width  between  scales  and  2.  The 

dyadic  scale  of  an  octave  is  the  last  voice  of  the  octave 
(scale  =  2°^^'"®^^).  In  Fig.  2,  we  consider  a  signal  of  512 
points  (2^).  This  means  9  octaves  (octave  0  to  octave 
8).  The  coarsest  scale  is  512,  the  finest  is 
1 

Y  _l_  ^  cxamplc,  when  we  display  a 

second  voice  of  the  fourth  octave  (four  voices  per 

2 

octave  computed),  we  obtain  the  scale  2^^^,  that  is  to 
say  scale  23. 

2.2.  One-dimensional  experiment 

We  applied  programs  from  libraries  in  LastWave  and 
Matlab,  using  a  continuous  wavelet  transform  and  a 
discrete  wavelet  transform.  LastWave  is  a  wavelet 


signal  and  image-processing  environment,  written  in  C 
[4].  Wavelab  is  an  extension  of  Matlab.  For  the  CWT, 
we  concentrated  on  the  first  and  the  second  derivative  of  a 
gaussian  function  (Mexican  Hat  wavelet).  We  processed 
phantom  signals  with  three  masses  of  distinct  sizes  using 
gaussian  additive  noise. 


Fig.  1.  The  time-frequency  lattice:  (a)  for  the  dyadic  wavelet 
transform,  is  localized  around  2”hbo;  ao=2  and  we  assume 
bo=l;  (b)  for  a  scheme  with  four  voices,  the  different  voice 
wavelets  are  assumed  to  be  dilatations  of  a  single 

function  ==  2^ ^^(2"^  x)  • 


section  of  a  mammograinm  without  mass 
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Fig.  4.  ID  sections  of  a  real  mammogram. 

As  shown  in  Fig.  4,  we  added  gaussian  noise  on  the 
phantom  mass  so  that  the  ID  signal  had  approximately 
the  same  shape  as  a  real  mass.  Fig.  5  shows  this 
representation. 
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Fig.  5.  Phantom  signal  with  an  added  white  gaussian  noise  of 
variance  0.1. 
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Fig.  2.  The  time-frequency  lattice  for  a  scheme  with  four  voices 
per  octave,  including  the  scale  axis. 

We  plotted  two  scan  line  profiles  (Fig.  4)  of  a  real 
mammogram  (Fig.  3). 
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Fig.  3.  Real  mass  from  a  mammogram.  The  white  lines 
show  the  locations  of  the  extracted  profiles  corresponding 
to  Fig.  4. 


Fig.  6  depicts  the  results  obtained  without 
downsampling.  The  signal  was  composed  of  masses 
with  a  white  gaussian  noise  of  variance  0.1.  The 
wavelet  was  a  Mexican  Hat. 
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Fig.  6.  Analysis  by  a  DWT  on  a  phantom  with  gaussian  noise 
(var  =  0.1).  We  show  approximation  and  detail  signals  at 
scales  8  (left)  and  16  (right). 


2.3.  The  2D  CWT 


We  began  the  2D  study  with  phantom  masses  of  white 
objects  on  a  black  background  with  the  addition  of  white 
gaussian  noise  of  variance  4.  We  applied  a  bias  to  the 
magnitude  values  to  preserve  the  waveform  shape  and 
make  the  signal  purely  positive  [5].  Next,  we  performed 
the  analysis  on  a  cancerous  mass  from  a  mammogram 
(Fig.  3).  We  show  the  biased  unthresholded  results  in 
Fig.  7,  and  the  thresholded  values  in  Fig.  8. 
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Fig.  7.  CWT_2D  at  octaves  4  to  6,  four  voices  per  octave.  No 
thresholding  on  biased  coefficients. 
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where  A“  (z)  is  the  autocorrelation  filter  of  degree  CC . 

The  transform  is  computed  for  a  real  mass 
along  the  scales  for  different  values  of  the  spline 
parameter  a  (Fig.  9). 
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Fig.  9-  DWT  in  2D  at  scale  4  and  8,  for  4  values  of  a. 


As  shown  in  Fig.  9,  we  do  not  always  observe  a  good 
representation  for  different  values  of  the  parameter  a. 
However,  we  clearly  observe  that  the  detection  is  better 
for  a=0.2.  The  parameter  of  the  spline  is  continuous 
(a>-0.5).  Therefore,  it  is  interesting  to  make  the 
parameter  vary  in  order  to  find  the  best  basis,  which 
suits  well  a  given  mass  size.  However  the  present 
transform  is  only  computed  at  dyadic  scales.  With  a 
continuous  analysis,  which  would  allow  decon:q)osition 
on  voices  between  these  scales,  we  may  obtain  a  richer 
parameter  space  so  as  to  identify  a  best  basis  for  mass 
detection. 


Fig.  8.  CWT_2D  at  octaves  4,  Sand  6,  four  voices  per  octave. 
Coefficients  are  biased  and  thresholded  among  scales  (10  for 
scale  38  to  20  for  scale  128). 

2.4.  Fractional  Splines 

We  have  more  recently  focused  on  the  Fractional  Spline 
Wavelet  Transform  [6,7].  We  have  extended  the 
implementation  to  two  dimensions,  which  was  described 
originally  by  M.  Unser  and  T.  Blu  [8].  We  used 
orthonormal  filters  to  compute  the  details  (horizontal, 
vertical  and  diagonal)  and  the  approximation  coefficients 
of  the  image  by  applying  the  filters, 


3.  Correlation  Analysis 

Given  the  results  in  one  dimension,  we  then 
implemented  a  2D  continuous  wavelet  transform.  Our 
goal  was  to  now  find  the  most  suitable  scale  to  detect  a 
mass  of  arbitrary  size.  To  find  the  best  scale,  we 
displayed  the  maxima  of  the  coefficients  along  scales, 
the  third  dimension  giving  the  magnitude  of  the 
maxima  at  each  scale.  In  addition,  we  plotted  the 
correlation  between  the  original  mass  and  the 
coefficients  of  the  CWT  at  each  scale.  We  expected  to 
find  different  optimal  scales  according  to  the  size  of  a 
mass.  We  tested  this  by  carrying  out  our  algorithm  on 
three  different  size  masses.  We  first  con:q)uted  for  each 
mass  the  CWT  in  2D  on  9  possible  octaves  (3  voices 
per  octave).  Then  for  each  octave  and  scale  we  plotted 


the  maxima  of  the  coefficients  of  the  wavelet 
decomposition  as  shown  in  Fig.  10. 


Fig,  10,  Evolution  of  the  maxima  of  the  cwt2d  across  scales. 

The  positions  of  the  maxima  of  the  decomposition  were  at 
scales  40,  81  and  128  for  small,  medium  and  large  masses 
respectively. 

Next,  we  performed  the  CWT  in  2D  on  the  same  number 
of  octaves  and  voices.  For  each  scale  we  calculated  the 
correlation  between  the  original  image  without  noise  and 
the  2D  CWT  decomposition  as  shown  in  Fig.  1 1 . 


Fig.  11.  Correlation  between  the  original  image  and  the  biased 
values  of  the  decomposition. 

The  positions  of  the  maxima  of  the  correlation  were  at 
scales  64,  102  and  161  for  small,  medium  and  large 
masses  respectively. 

The  most  suitable  scale  using  the  method  of  the  maxima 
evolution  was  not  the  same  as  the  scale  identified  with 
correlation.  Next,  we  attempted  to  find  the  best  scale  for 
a  real  mass.  This  time,  the  best  scale  to  detect  the  real 
mass  was  the  same  for  both  methods  (maxima  evolution 
and  correlation)  at  scale  161.  We  also  considered  a  very 
noisy  signal  (variance  4),  for  robustness.  We  analyzed 
the  maxima  of  the  coefficients  and  the  correlation  for 
different  noise  settings.  From  these  results  we  observed 
that  both  methods  identified  same  scale  values  regardless 
of  the  amoimt  of  added  noise. 


4.  Multi-scale  Adaptive  Histogram  Equalization 

Analog  and  digital  mammography  often  contains  12 
bits  or  more  of  significant  contrast  information. 
Anatomical  tissues  may  occupy  significantly  different 
dynamic  ranges  on  display  due  to  difference  of  X-ray 
attenuation.  By  comparison,  the  human  visual  system 
can  only  perceive  less  than  100  different  gray  levels 
[9].  Thus,  contrast  enhancement  is  usually  needed  for 
clinical  readings.  This  section  discusses  one  approach 
for  contrast  enhancement  utilizing  multi-scale  analysis. 
Sub-band  coefficients  were  modified  by  the  method  of 
adaptive  histogram  equalization.  To  achieve  optimal 
contrast  enhancement,  the  sizes  of  sub-regions  were 
chosen  with  consideration  to  the  support  of  the  analysis 
filters.  The  enhanced  images  provided  subtle  details  of 
tissues  that  are  only  visible  with  tedious 
contrast/brightness  windowing  methods  currently  used 
in  clinical  reading. 

By  properly  selecting  the  decomposition  filters,  desired 
features  of  an  object  can  be  separated  from  noise. 
Therefore  we  can  selectively  enhance  features  of 
interest  by  modifying  corresponding  components  in  the 
transform  domain.  We  used  the  quadratic  spline 
wavelet  function  y/{x) ,  which  has  compact  support 
and  is  continuously  differentiable.  It  is  the  derivative 
of  a  cubic  spline  function  6{x)  as  seen  in  Fig.  12.  It 
can  be  shown  that  by  using  a  wavelet  that  is  the 
derivative  of  a  smoothing  function  the  wavelet 

transform  f  of  the  signal  f  is  proportional  to  the 

derivative  of  the  signal  smoothed  at  scale  2*^ .  The 
wavelet  transform  can  then  be  considered  as  an 
adaptive  (scale  dependent)  detection  procedure  that 
finds  signal  variation  points  in  two  orthogonal 
directions  x  andy  [10]. 


(a)  (b) 


Fig  12.  (a)  Cubic  spline  smoothing  fiinction  6(x).  (b) 
Quadratic  spline  wavelet  xf/(x)  of  compact  support  defined  as 
the  derivative  of  the  smoothing  function. 

In  Fig  13,  we  present  results  on  mammography  data, 
which  shows  significant  improvement  over  existing 
traditional  window  and  leveling  techniques  used  in 
soft-copy  stations.  The  contrast  limited  adaptive 
histogram  equalization  (CLAHE)  clearly  enhances 
monographic  features. 


(b)  (c) 


Fig.  13.  (a)  original  film,  (b)  detailed  window  & 
leveling  by  a  radiologist,  (c)  contrast  limited  adaptive 
histogram  equalization  (CLAHE). 


5,  ROC  Study:  Experimental  Design 
Our  study  focused  on  density  3  and  4  mammograms,  on 
BiRads  scale,  A  total  of  60  cases,  subdivided  into  2 
groups  (15  cancers  &  15  normals)  were  read  by  three 
radiologist  [11],  The  diagnosis  for  mammograms 
included  BI-RAD  (0-5),  LOG  value  (1-5)  and 
localization  of  detected  lesion. 

We  have  divided  the  radiologist  into  2  groups: 

Group  1:  “Softcopy  Display”  +  Interactive  Contrast 
Enhancement. 

Group  2:  “Softcopy  Display”  only. 

Our  softcopy  monitors  consist  of  high-resolution 
(2048x2560)  dual  Barco  5MP1H  displays.  Fig.  15. 
depicts  a  radiologist  participant  in  the  study  along  with 
dual  monitor  displaying  monographic  data. 


Fig.  15.  Actual,  experimental  study. 


The  ROC  analysis  can  be  seen  in  Fig.  16.  The  area 
under  the  curve  with  computer  aided  enhancement  and 
without  enhancement  is  0.9136  and  0,8405 
respectively.  The  computer-aided  diagnosis  brings  a 
noticeable  improvement  in  cancer  screening. 
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Fig.  16.  ROC  curves  for  data  from  Group  1  (“with 
Enhancement”)  and  from  Group  2  (“without  Enhancement”) 
[11]. 


6.  CONCLUSION 

Our  studies  in  one  and  two  dimensions  suggest  that 
dyadic  scales  are  often  not  sufficient  to  detect  a  mass  in  a 
dense  mammogram.  We  showed  the  advantage  of  a 
continuous  wavelet  transform,  which  computed  an 
expansion  on  voices  between  the  common  dyadic  scales. 
We  saw  on  real  images  of  masses  extracted  from  digitized 
mammograms  that  a  correlation  method  between  a  known 
mass  and  the  values  of  computed  coefficients  yielded 
approximately  the  same  results,  as  a  maximum  method 
evolution.  Thus,  this  study  suggests  that  it  is  possible  and 
of  value  to  tune  an  analysis  between  octaves,  for  the 
detection  of  subtle  masses  in  mammograms. 

On  a  second  front,  a  multi-scale  adaptive  histogram 
equalization  method  was  reported  here,  which  showed 
promising  results  on  mammography  interpretation.  We 
claimed  that  the  advantage  of  this  method  comes  from 
combining  the  local  enhancement  ability  of  AHE,  and  the 
selectivity  of  spatial-frequency  conqjonents  from  wavelet 
analysis.  The  overall  diagnostic  sensitivity  compares 
favorable  with  state-of-art  enhancement  methods,  and 
also  circumvents  and  reduces  some  of  the  artifacts 
visualized  with  existing  methods.  The  ability  of 
simultaneously  displaying  the  full  dynamic  contrast  range 
was  shown  to  be  efficient  in  terms  of  interpretation  time. 
The  diagnostic  performance  showed  the  possibility  of 
building  a  new  “Power  Windows”  scheme  for  clinical 
usage. 

Following  the  conventional  three-windows  settings,  we 
can  also  tailor  the  parameters  to  find  the  best 
enhancement  for  particular  abnormalities  based  on  their 
spatial-frequency  properties.  This  method  proves  to  be  of 
value  in  isolating  cancer  masses  of  diameter  1cm  or  less. 
We  certainly  expect  a  more  reliable  diagnosis  compared 
to  existing  windowing  schemes. 
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ABSTRACT 

We  have  carried  out  a  receiver  operating  characteristics  (ROC)  study  for  the  enhancement  of  mammographic  features  in 
digitized  mammograms.  The  study  evaluated  the  benefits  of  multi-scale  enhancement  methods  in  terms  of  diagnostic 
performance  of  radiologists.  The  enhancement  protocol  relied  on  multi-scale  expansions  and  non-linear  enhancement 
functions.  Dyadic  spline  wavelet  functions  (first  derivative  of  a  cubic  spline)  were  used  together  with  a  sigmoidal  non-linear 
enhancement  function  [1],  [2].  We  designed  a  computer  interface  on  a  softcopy  display  and  performed  an  ROC  study  with 
three  radiologists,  who  specialized  in  mammography.  Clinical  cases  were  obtained  from  a  national  mammography  database 
of  digitized  radiographs  prepared  by  the  University  of  South  Florida  (USF)  and  Harvard  Medical  School. 

Our  study  focused  on  dense  mammograms,  i.e.  mammograms  of  density  3  and  4  on  the  American  College  of  Radiology 
(ACR)  breast  density  rating,  which  are  the  most  difficult  cases  in  screening,  were  selected.  To  compare  the  performance  of 
radiologists  with  and  without  using  multi-scale  enhancement,  two  groups  of  30  cases  each  were  diagnosed.  Each  group 
contained  15  cases  of  cancerous  and  15  cases  of  normal  mammograms.  Conventional  ROC  analysis  was  applied,  and  the 
resulting  ROC  curves  indicated  improved  diagnostic  performance  when  radiologists  used  multi-scale  non-linear 
enhancement. 

Keywords:  Multi-scale  analysis,  ROC  analysis,  contrast  enhancement,  digital  mammography,  softcopy  display. 

1.  INTRODUCTION 

Recently,  research  has  focused  on  the  development  of  digital  displays  and  softcopy  workstations  for  digital  mammography. 
Limited  spatial  resolution,  luminance,  and  dynamic  range  cannot  be  solved  simply  by  hardware  improvements  or  computer 
programming  alone.  A  possible  solution  of  these  problems  is  the  application  of  multi-scale  contrast  enhancement  techniques 
derived  from  non-linear  models. 

Radiologists  are  mostly  familiar  with  films  where  the  Modulation  Transfer  Function  (MTF)  is  approximately  equal  to  2®  gray 
levels  of  contrast  resolution.  However,  images  acquired  with  digital  detectors  can  record  at  least  2*^  different  gray  levels  of 
intensity  and  are  now  commercially  available.  The  wealth  of  dynamic  range  within  these  digital  acquisition  systems  provides 
strong  evidence  that  the  signal-to-noise-ratio  (SNR)  can  be  increased  in  digital  mammography.  For  expert  radiologists  the 
human  visual  system  can  detect  at  most  2’  shades  of  gray.  These  considerations  motivate  the  need  for  judicious  methods  of 
processing  of  digital  radiographs  that  can  optimize  the  bandwidth  of  the  human  visual  system.  We  have  designed 
enhancement  software  that  is  well  adapted  for  this  purpose  and  provides  a  “data  mining”  tool  to  map  and  make  visible 
selected  “quantum  levels”  of  information  living  within  the  wide  range  of  contrast  resolution  provided  by  digital  detectors. 

Medical  imaging  is  a  field  in  which  quantitative  accuracy  and  qualitative  fidelity  are  paramount.  In  any  image  enhancement 
process  distortion  of  the  original  image  and  artifacts  are  not  affordable.  Multidimensional  feature  enhancement  via  wavelet 
analysis  has  been  previously  demonstrated  on  mammograms  [3],  [4],  [5],  [6],  [7],  [8]  and  is  a  powerful  tool  for  processing 
digital  medical  images  without  artifacts.  The  enhancement  process  adjusts  multi-scale  coefficients  at  some  particular  spatial- 
frequency  scale  by  increasing,  decreasing  or  resetting  their  values.  Each  image  is  then  reconstructed  with  modified 
coefficients.  This  simple  enhancement  technique  relies  on  the  idea  that  features  of  interest  in  a  given  radiograph  are 
detectable  at  a  particular  scale  and  can  be  amplified,  whereas  noise  and  less  clinically  interesting  features  may  live  at  other 
levels  of  analysis  whose  visual  appearance  can  be  diminished  or  eliminated  in  a  reconstructed  image.  Further  results  and 
detailed  descriptions  of  these  methods  can  be  found  in  [9],  [10],  [11],  [12],  [13],  [14],  [15]. 

Surprisingly,  there  have  been  very  few  studies  carried  out  to  evaluate  the  benefits  of  multi-scale  enhancement  methods  in 
terms  of  diagnostic  performance.  Our  study  aimed  at  providing  quantitative  evidence  of  these  benefits.  ROC  analysis  [16]  is 
most  commonly  used  in  medical  imaging  for  such  purposes,  though  alternative  statistical  approaches  can  be  found  as  well 
[17].  ROC  curves  have  been  compared  to  evaluate  the  visibility  of  malignancies  [18],  mass  detection  techniques  [19]  or 
algorithms  for  computer-aided  diagnosis  (CAD)  that  use  neural  networks  [20], 


The  chapter  is  organi^d  as  follows.  In  Section  2  we  describe  a  protocol  for  multi-scale  non-linem*  contet  enhancement. 
After  a  short  overview  of  the  use  of  multi-scale  expansions  for  contrast  enhancement  we  discuss  the  dyadic  spline  wavelet 
selected,  its  implementation,  and  how  a  non-linear  enhancement  function  is  applied  to  multi-scale  coefficients.  Section  3 
addresses  the  design  of  a  graphical  user  interface  (GUI)  that  was  developed  to  caixy  out  the  ROC  study  including  high- 
performance  displays  and  specialized  hardware  for  softeopy  display  of  digital  mammograms.  Next,  ftie  ROC  study  itself 
together  with  iU  lesults  and  subsequent  data  analysis  is  pmsented  in  Section  4.  After  a  discussion  of  the  resulte  of  the  study, 
conclusions  and  possible  directions  of  fliture  research  are  presented  in  Section  5. 


2.  ENHANCEMENT  PROTOCOL 


2.1.  Contrast  Enhancement  via  Multi-scale  Expansions:  A  Short  Overview 


We  summarize  below,  the  advantages  of  the  use  of  overcomplete  multi-scale  representotions  for  adaptive  contrast 
enhancement  of  digital  mammograms.  Critically  sampled  multi-scale  representations  are  not  suitable  for  detection  and 
enhancement  tasks  because  of  alining  effecte  introduced  during  downsampling  of  the  analysis  [21],  [22].  However, 
overcomplete  representations  avoid  such  aliasing  artifacts  and  have  the  desirable  properly  of  being  shift  invariant  [23],  [24]. 
Indeed,  this  property  ensures  that  die  spatial  locations  of  any  mammographic  finding  wifliin  in  an  image  are  preserved  across 
all  scales.  Thus,  in  our  approach  the  transform  coefBcient  matrix  si^  at  each  scale  remains  the  same  as  the  original  spatial 
resolution  of  die  digital  mammogram,  since  there  is  no  downsampling  across  each  level  of  analysis. 


Overcomplete  multi-scale  analysis  and  leconsfruction  algoritfims  using  dyadic  scales  previously  developed  in  [25],  [26],  md 
[27]  were  used  as  an  initial  choice  of  analysis  function  for  our  enhancement  protocol.  The  implementetion  was  cmried  out 
using  several  lowp^s  and  highpass  filters  with  localized  frequency  support.  At  each  level  of  die  multi-scale  expansion  m 
input  image  is  decomposed  into  a  coarse  approximation  and  detailed  structures.  The  coarse  approximation  is  the  ou^ut  from 
applying  a  lowpass  filter,  and  the  detailed  structures  are  obtained  fi-om  highp^s  filtering.  The  approximation  image 
corresponds  to  scaling  coefficients,  whercas  the  details  extracted  fiom  the  approximation  are  wavelet  coefficiente  at  a 
particular  scale.  This  procedure  is  successively  repeated  on  the  ^proximation  image  to  obtain  multiple  levels  of  analysis, 
Ihe  coaraest  approximation  is  often  referred  to  as  “dc-cap”.  A  gain  or  enhancement  ftmction  modifies  the  matrices  of 
coefficiente  that  have  been  isolated  by  the  filtera  at  each  level  and  may  boost  coefficients  at  some  scales  and/or  attenuate 
others.  If  the  filters  meet  a  perfect  reconstniction  condition,  toe  image  can  be  reconstructed  from  its  wavelet  representation  of 
scaling  and  wavelet  coefficients  [28].  The  filter  bank  implementation  of  enhancement  processing  by  an  expansion- 
reconstiuction  algorithm  for  2  levels  of  analysis  is  schematically  illustrated  in  Figure  1,  Image  rccx)nstruction  that  is  also 
accomplished  by  appropriate  filtering  operations  is  presented  in  a  simplified  manner  in  Figure  1. 
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Figure  1:  Multi-scale  analysis  with  non-linear  contrast  enhancement:  Schematic  of  filter  bank  implementation.  In  the  left  part 
multi-scale  expansion  with  enhancement  for  2  levels  of  analysis  is  shown,  and  raconstruction  is  presented  (in  a  simplified 
manner)  in  the  right  part 

The  modified  matrices  of  coefficients  are  simply  ‘‘plugged  in”  during  reconstruction  producing  a  “focused”  subband 
enhancement.  As  shown  above,  the  enhancement  fimction  can  be  implemented  independently  of  a  particular  set  of  filtera  mid 
easily  incorporated  into  a  filter  bank  to  provide  the  benefite  of  multi-scale  enhancement  [1],  [29]. 


2.2.  High  Speed  Implementation  to  Support  Interactive  Processing 

Similar  to  orthogonal  and  biorthogonal  discrete  wavelet  transforms  [30],  the  discrete  dyadic  wavelet  transform  can  be 
implemented  within  a  hierarchical  filtering  scheme.  Let  an  input  signal  x(n)  be  real,  x(n)  g  l\Z),  «  e  [O,  i.e.,  x(n)  is 

supported  on  the  index  interval  [0,  N-1],  and  let  X{g))  be  its  Fourier  transform.  Depending  on  the  length  of  each  filter 
impulse  response,  filtering  an  input  signal  may  be  computed  either  by  multiplying  X(o))  by  the  frequency  response  of  a 
filter  or  by  circularly  convolving  x(n)  with  the  impulse  response  of  a  filter.  Of  course,  such  a  periodically  extended  signal 
may  change  abruptly  at  the  boundaries  and  cause  artifacts.  A  common  remedy  for  such  a  problem  is  realized  by  constructing 
a  mirror  extended  signal 


x{n). 


where  we  chose  the  signal  Xj„e(n)  to  be  supported  in  [-N,  N-1].  In  [1]  it  is  shown  how  a  mirror  extension  is  a  particularly 
elegant  solution  in  conjunction  with  symmetric/anti-symmetric  filters,  since  a  signal  is  of  a  particular  type  of  symmetry  at 
each  stage  of  the  filter  bank.  The  optimized  circular  convolution  described  in  [1]  was  implemented  in  native  “ANSI  C”  to 
speed  up  performance  for  multi-scale  decomposition  and  image  reconstruction.  Parameters  of  this  algorithm  included  number 
of  levels  of  analysis,  gain,  and  threshold.  This  algorithm  was  incorporated  into  a  graphical  user  interface  (GUI)  developed 
during  the  preparation  of  the  study. 


As  a  further  goal,  we  envision  developing  feature  specific  enhancement  protocols  for  each  type  of  lesion.  An  enhancement 
protocol  would  consist  of  a  multi-scale  expansion  of  a  mammogram  by  a  specific  basis  and  an  associated  non-linear 
enhancement  function  that  is  best  matched  to  a  specific  type  of  lesion,  e.g.  microcalcifications.  For  the  study  under 
consideration,  a  dyadic  spline  wavelet  function  was  used  as  the  basis,  and  a  non-linear  sigmoidal  function  was  applied  as  the 
enhancement  fimction.  Both  are  described  in  greater  detail  next. 

2.3.  Dyadic  Spline  Wavelet  Algorithm 


The  wavelet  transform  of  a  signal  f(x)  at  scale  5*  and  position  x  is  defined  by  Wf(u,s)  =  /  *  ^  “  JC  ^  ’ 

where  the  function /is  projected  on  a  family  of  translated  and  dilated  basis  functions  (wavelets)  ^„,(x)  =  i^[  ^  1. 

K  s  ) 

y/{x)  is  the  mother  wavelet  of  zero  average.  Both,  translation  and  dilation  parameters  u  and  s  are  continuous  for  the 
continuous  wavelet  transform.  To  allow  fast  numerical  implementation  of  discrete  wavelet  transforms,  Mallat  and  Zhong  [31] 
introduced  the  dyadic  wavelet  transform,  where  the  scale  parameter  varies  only  along  the  dyadic  sequence  {2^},  with  j  eZ . 

Extending  this  approach  to  two  dimensions  by  the  use  of  a  tensor  product  yields  the  2-D  dyadic  wavelet  transform  that 
partitions  plane  orientations  into  two  bands.  This  means  that  there  are  two  channels  of  analysis  along  orthogonal  directions  x 
and  y.  The  wavelet  transform  of  the  2-D  signal  f(x,y)  at  the  scale  2^  has  two  components  defined  by: 

Wljf{x,y)  =  f*ii/\j{x,y)  and  W^,f{x,y)  =  f*y/l,{x,y),yi\ih.  =  (d=I,2).  We  used  the  quadratic 

spline  wavelet  function  ^(x)  defined  by  Mallat  and  Zhong  in  [31]  of  compact  support  and  continuously  differentiable.  Its 


Fourier  transform  can  be  derived  as  y/{p))  =  {jm) 

function  0(x) ,  whose  Fourier  transform  is  S(co)  = 
case  in  Figure  2  below. 
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i//{x)  is  the  first  derivative  of  a  cubic  spline  smoothing 
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[1].  These  functions  are  displayed  for  the  one-dimensional 


(«)  (b) 

Figure  2t  (a)  Cubic  spline  smoothing  function  $(x).  (b)  Quadratic  spline  wavelet  ^x)  of  compact  support  defined  as  the  ilist 
derivative  of  the  smoothing  function. 

Using  a  wavelet  tiiat  is  the  derivative  of  a  smoothing  fimction  it  can  be  shown  that  the  wavelet  transform  of  the  signal / 

is  proportional  to  tiie  derivative  of  the  signal  smoothed  at  the  scale  2^  [32].  The  coefficients  of  modulus  maxima  detection  are 
then  equivalent  to  an  Captive  sampling  that  finds  signal  variation  pointe  in  the  two  orthogonal  directions  x  and  y. 

As  images  represent  finite  enei^  signals  measured  at  some  finite  ^solution,  we  cannot  compute  the  wavelet  transform  at 
scales  below  the  limit  set  by  this  resolution.  We  applied  this  analysis  at  dyadic  scales  varying  fiom  1  (original  signal)  to  the 
limit  imposed  by  acquisition  (digiti^r  sampling  rate).  Figure  3  shows  an  example  for  one  level  of  an  overcomplete  wavelet 
expansion  of  a  region  of  interest  (ROI)  with  a  spiculated  m^s  at  a  dyadic  scale,  and  in  Figure  4  wavelet  coefficients  of 
microcalcifications  at  the  finest  dyadic  scale  are  presented. 


(a)  (b)  (c)  (d) 


Figure  3:  Level  5  of  an  overcomplete  dyadic  wavelet  expansion  of  a  spiculated  mass,  (a)  Original  image.  Q>)  Horizontal  details, 
(c)  Vertical  detail$.(d)  Approximation  image. 


(a)  (b)  (c) 

Figure  4:  (a)  Original  ROI  with  mkrocalcifications.  (b)  Horizontal  and  (c)  Vertical  dyadic  wavelet  coefficients. 

2.4*  Non-Linear  Enhancement  Function 

Modification  of  selected  analysis  coefficients  within  a  certain  scale  can  make  more  obvious  indiscernible  or  barely  seen 
mammographic  features  [14].  Contrast  enhancement  was  achieved  by  applying  an  enhancement  fimction  to  transform 
coefficients  at  selected  scales.  This  operation  results  in  local  attenuation  or  amplification  of  coefficients.  Enhancement  or 
gain  fimctions  must  be  cumulative  and  monotonically  increasing,  in  order  to  preserve  the  order  of  intensity  information  in  the 
original  image  and  to  avoid  artifacts  [26].  Figure  5(a)  provides  a  very  simple  example  of  a  piecewise  linear  enhancement 
function.  Multi-scale  coefficients  are  denoted  Wy,  which  are  modified  by  applying  an  enhancement  fimction T  is  tiie 
threshold  of  tiie  function,  and  a  foe  gain.  The  effect  of  foe  enhancement  function  depends  on  the  value  of  the  angle  0,  For  0  < 
45*  there  is  an  attenuation  of  foe  coefficients  (a<l),  at  ^  =  45*  we  have  foe  identity  function  (a=l),  and  for  0  >  45*  foere  is  a 
smooth  amplification  of  foe  coefficients  (^1).  The  values  of  foe  two  parameters,  T  and  0  (or  a),  determine  the  final  shape  of 
foe  enhancement  fimction.  Figure  5(b)  displays  a  hard-thresholding  function  for  denoising,  where  coefficiente  wifo  modulus 
<  r  are  set  to  zero.  Unfortunately,  these  two  particular  fimctions  have  the  disadvantage  of  being  discontinuous  at  foe 
threshold  value  :ff.  This  could  result  m  an  abnormal  distribution  of  coefficient  values  in  foe  output  and  may  create  sharp 


peaks  on  both  ends  of  the  histogram  of  a  particular  output  mapping.  For  this  reason,  smoother  functions,  like  sigmoids,  are 
preferable  and  were  used  in  this  study.  Figure  6  shows  an  example  of  such  a  function  as  described  in  [2]. 


Figure  6:  A  sigmoidal  non-linear  enhancement  function. 

The  analytical  formulation  of  the  sigmoidal  enhancement  function  as  designed  in  [33],  [2]  is  the  following: 

/( Wjj )  =  a  [sigm  (c(Wij  -b))-  sigm  (-c(Wij  +  i))] 

a  = - z - r - - - r,  0<Z)<1 

sigm  (c(l  -  Z)))  -  sigm  {-c(l  +  Z))) 

sigm{y)  =  —^ 

\  +  e  ^ 


(1) 


Parameters  b  and  c  control  the  threshold  and  the  rate  of  enhancement  (gain)  respectively.  This  enhancement  function  is 
continuous,  monotically  increasing,  and  has  a  continuous  first  derivative.  This  ensures  that  the  application  of  the  function  will 
not  introduce  any  new  discontinuities  of  coefficients  in  the  transform  domain. 

From  Figure  6  we  see  that  this  enhancement  function  decreases  the  value  of  the  coefficients  around  zero,  which  is  equivalent 
to  a  denoising  action,  while  it  may  increase  values  of  the  coefficients  outside  this  range,  equivalent  to  enhancement  or 
amplification.  This  type  of  enhancement  function,  in  ‘steps’,  offers  a  very  rich  and  flexible  paradigm  to  carry  out  non-linear 
dynamic  analysis  of  coefficients  within  a  specific  scale  [34]. 

There  are  many  criteria  for  the  selection  of  the  enhancement  function  applied  to  the  coefficients  of  a  particular  level  of 
analysis  for  contrast  enhancement.  One  goal  of  the  study  described  here  was  to  develop  a  research  tool  for  testing 
enhancement  functions  targeted  for  specific  mammographic  features.  As  this  process  requires  specialized  expertise  and  a 
substantial  time  investment,  no  systematic  study  of  the  problem  of  associating  enhancement  functions  with  target  features  in 
mammograms  has  been  reported  in  the  literature. 

In  general,  non-linear  estimators  are  signal  dependent  and  behave  differently  for  different  realizations  of  each  signal.  In  this 
context,  Johnstone  and  Donoho  have  shown  that  by  considering  the  signal  as  deterministic,  thresholding  of  wavelet 
coefficients  gives  a  nearly  optimal  estimation  of  piecewise  smooth  functions  [35],  [36].  More  specifically,  for  a  noisy  signal 

of  size  N,  thresholding  of  the  wavelet  coefficients  with  T  =  (7^j2\n(N) ,  where  a  is  the  standard  deviation  of  the  coefficients, 
provides  an  asymptotically  optimal  estimator  of  the  original  signal  in  the  mini-max  sense  [36].  Thresholding  of  wavelet 
coefficients  performs  an  adaptive  smoothing  of  the  image  by  averaging  noisy  areas  and  preserving  or  enhancing  coefficients 
in  areas  of  sharp  transitions.  Noise  standard  deviations  can  be  estimated  by  determining  the  median  wavelet  coefficient  value 
at  the  finest  scale  or  with  local  discrete  statistical  estimation  in  the  transform  domain.  Using  extremely  local  variances  for  the 


estimation  of  a  tiireshold  leads  to  a  very  aggressive  posturing  of  the  enhancement  functions  snd  represents  a  high  amount  of 
intervention  in  mijusting  the  output,  while  global  variance  measurements  are  less  noticeable.  Superiority  of  either  method 
depends  on  the  screening  protocol  used  by  the  radiologist  and  the  kind  of  analysis  to  be  performed.  For  example,  fine 
microcalcifications  represent  high  frequency  information  of  foe  image.  We  would  expect  foe  local  variance  for  such  a  feature 
to  be  high  within  a  selected  ROL  Consequently,  smooth  amplification  of  coefficients  within  this  particular  spatial  frequency 
range  (in  combination  with  possibly  decreasing  foe  information  of  ofoer  spatial  frequencies)  will  enhance  foese  features  of 
interest.  Similar  analysis  can  be  done  to  enhance  low  spatial  frequency  features  such  as  masses.  A  block  diagram  of  the 
enhancement  process  for  coefficients  at  selected  scales,  which  me  chosen  with  respect  to  the  particular  mammographic 
feature  to  be  danced,  is  shown  in  Figure  7  below. 

Since  the  computation  of  the  enhancement  parameters  uses  data  dependent  information  such  as  local  or  global  coefficient 
variance,  digital  and  digitized  radiographs  acquired  imder  different  imaging  conditions  are  best  processed  independently  to 
achieve  optimal  enhancement.  Intrinsic  properties  of  foe  radiograph  are  therefore  incorporated  in  foe  setting  of  foe 
paimieters.  In  our  work  we  used  both  coefficient  variance  compute  with  respect  to  a  selected  ROI  mid  user  input  (see 
Section  3.2)  to  adapt  the  foreshold  and  gain  parameters. 


Figure  7:  Block  diagram  of  modifying  feature  specific  coefficients  at  selected  scales  by  appb^lng  a  non-linear  enhancement  function. 

3.  DEVELOPMENT  OF  A  GRAPHICAL  USER  INTERFACE  (GUI) 

3.1.  Motivation 

Ruiming  such  an  erfoancement  algorithm  in  a  batch  mode  might  be  sufficient  for  single  experiments.  However,  adjustment  of 
parametera  tied  to  a  data  dependent  enhancement  fimetion  is  slow,  because  of  the  repeated  need  to  decompose  ^d 
reconstruct  from  modified  coefficients.  A  more  desirable  situation  would  be  to  observe  foe  resulte  of  modified  multi-scale 
croefficiente  interactively  mid  to  crontinue  the  enhancement  procedure,  until  results  are  visually  satisfactory  or  the  decision  is 
made  that  no  further  improvement  can  be  achieved.  In  addition,  with  introducing  fixed  enhancement  protocols  into  a  clinical 
screening  paradigm,  the  algorithm  must  be  simple,  fast,  and  user-friendly,  i.e.  usage  of  foe  algorithm  should  be  femiliar  to  foe 
radiologist  and  intuitive.  Since  each  radiologist  may  have  preferences  with  respect  to  mnUmt  in  mammograms,  it  must  be 
IKissible  to  mijust  parameter  settings  to  individual  preferences.  Thus,  we  designed  a  graphical  user  interface  (GUI)  to 
facilitate  carrying  out  such  a  study  and  to  create  a  softcopy  display  prototype,  whose  successors  might  find  entrance  into 
clinical  screening.  We  call  this  application  a  ‘*test  bed”  softcopy  display  tool.  Its  first  veraion  was  employed  for  foe  ROC 
study  described  in  foe  next  s^tion. 

3.2.  Design  and  Implementation 

The  graphical  user  interface  (GUI)  developed  for  this  study  w^  written  in  Visual  C++  6.0.  The  code  for  foe  wavelet 
expansion  and  image  reconstruction  that  was  written  in  native  *‘ANSI  C”  to  speed  up  performance  could  be  incorporated  and 
executed  in  this  environment  without  major  modifications,  thus  shortening  development  time.  Some  of  the  guidelines  and 
considerations  for  foe  design  and  implementation  of  foe  GUI  are  described  next. 

The  prototype  interface  was  primarily  designed  to  process  raw  16-bit  data.  Data  was  obtained  from  a  national  mammography 
database  of  digiti^d  radiographs  provided  by  the  University  of  South  Florida  (USF,  ‘^Digital  Database  for  Screening 
Mammography”  (DDSM)).  Our  database  of  digitized  mammograms  (stored  on  twenty-two  8mm  tapes)  at  foe  time  of  foe 
study  contain^  586  selected  cases  of  malignant  lesions,  biopsy  proven,  and  437  cases  of  normal  breasts.  More  specifically, 
different  types  of  lesions  are  represented  in  foe  following  proportions:  100  round  and  oval  malignant  masses,  216  spicular 
lesions  and  248  microcalcifications.  559  c^es  of  dense  breaste  (density  of  3  and  4)  with  266  normals  and  293  cancerous, 
referred  by  radiologists  as  the  most  challenging  cases,  were  included  in  the  database. 

Images  from  foe  mammography  database  were  digiti^d  from  film  at  resolutions  of  40  to  50  pm.  Im^e  line  lengths  (#  of 
columns)  varied  between  2000  and  3000  pixels,  and  number  of  rows  from  4000  to  5900  pixels.  Depending  on  the  scanner 
utilized  for  digitization  foe  contrast  resolution  was  either  12  bits  or  16  bits  per  pixel  resulting  in  15-50  megabytes  per  view. 

To  handle  this  large  amount  of  data  and  to  provide  foe  diagnosing  radiologist  as  much  information  as  possible,  all  four  views 
(right  and  left  medial-lateral  (RMLO,  LMLO)  and  right  and  left  cranial-caudal  (RCC,  LCC)  of  a  case  were  loaded  into 
memory  and  displayed  as  downsampled  images  on  display  screen,  which  consisted  of  two  high-resolution  MegaScan 


monitors  each  with  a  screen  size  of  2048  by  2560  pixels.  Specialized  framebuffers  allowed  a  display  of  2^°  gray  levels  (see 
Section  3.3).  The  four  views  were  aligned  to  assist  the  radiologist  to  look  for  asymmetries.  In  addition,  one  view  could  be 
selected,  and  a  viewport  could  display  a  selected  region  of  interest  (ROI)  at  full  (original)  resolution  from  a  selected 
mammogram.  The  size  of  the  viewport  could  be  chosen  as  512  by  512,  1014  by  1024  or  even  2048  by  2048.  The  center  of  the 
ROI  was  determined  through  the  mouse  pointer  in  a  chosen  window.  Thus,  the  original  mammogram  could  also  be  examined 
through  the  viewport,  if  desired.  More  importantly,  suspicious  areas  could  be  captured  in  the  viewport  and  processed  through 
enhancement  via  the  multi-scale  expansion  described  in  Section  2.  For  the  enhancement  procedure  the  user  could  adjust  the 
number  of  subbands  of  the  expansion  as  well.  After  selecting  a  ROI  the  image  was  decomposed  onto  dyadic  wavelet  basis 
functions  yielding  wavelet  coefficients.  Coefficients  were  modified  by  a  sigmoidal  non-linear  enhancement  function,  and  the 
image  was  reconstructed  from  these  modified  coefficients  in  nearly  real-time. 

Figure  8(a)  shows  Dr.  Koenigsberg,  one  of  three  radiologists  who  participated  in  this  investigation,  during  the  ROC  study. 
Figure  8(b)  depicts  a  typical  screen  display  of  the  GUI  showing  additional  viewports  described  above. 


(a)  (b) 


Figure  8:  (a)  Tova  Koenigsberg,  M.D.,  using  the  GUI  during  the  preliminary  ROC  study  described  above,  (b)  Typical  screen 
display  used  during  the  ROC  study:  four  original  digitized  mammograms  of  one  case  on  the  right  monitor,  and  a 
selected  view,  the  GUI  interface  for  parameter  adjustments,  original  and  enhanced  ROI  are  shown  on  the  left  monitor. 

As  mentioned  in  Section  2.4  the  shape  of  the  enhancement  function  can  be  changed  through  modification  of  the  two 
parameters  gain  and  threshold.  Therefore,  each  parameter  could  be  adjusted  through  sliders  for  each  level  (subband)  of  the 
multi-scale  expansion  (see  Figure  9(b)).  On  release  of  the  slider  button,  a  reconstruction  “event”  was  “triggered”,  and  a 
resulting  image  presented  in  an  output  window.  For  example,  reconstruction  ofa512by512  matrix  for  five  levels  of 
decomposition  (5  subbands)  took  5  to  6  seconds.  For  four  subbands  reconstruction  time  shortened  to  4  to  5  seconds. 
Reconstruction  times  trecon  for  different  sizes  of  the  ROI  and  different  number  of  levels  of  analysis  are  presented  in  Table  1. 
However,  reconstruction  time  can  certainly  be  improved  to  achieve  true  real-time  performance,  by  employing  faster 
algorithms. 


512x512 

4-5  seconds 

6-7  seconds 

1024  X 1024 

19-20  seconds 

24-25  seconds 

Table  1:  Reconstruction  times  trecon  for  two  different  levels  of  analysis  and  two  sizes  of  ROI. 


After  processing,  enhanced  images  could  be  saved  together  with  information  about  the  location  of  the  ROI  (the  position  of 
the  ROI  was  marked  in  its  corresponding  downsampled  view)  to  facilitate  evaluation  of  a  particular  diagnosis  for  each  case  in 
comparison  with  the  “ground  truth”  provided  in  the  USF  database.  All  suspicious  areas  in  a  case  could  be  carefully  examined 
by  sequentially  choosing  different  views  and  multiple  ROIs. 

Figure  9(b)  shows  the  test  bed  interface  as  an  illustration.  Interactive  (real-time)  enhancement  was  accomplished  via  sliders 
shown  in  the  graphical  user  interface  (GUI).  The  enhancement  operation  relied  on  the  optimality  of  parameters  derived  from 
their  non-linear  models  and  on  the  strategy  employed  for  the  type  of  enhancement  applied  to  each  subband  of  coefficients 
(amplification,  preservation  or  diminution).  Selected  subband  coefficients  at  a  particular  level  could  be  strongly  suppressed 
by  choosing  large  thresholds  (>  2)  and  small  gains  (<  1),  which  can  be  desirable  for  the  elimination  of  (structured  and 
acquisition)  noise,  or  normal  benign  anatomical  (fibroglandular)  structures. 

Since  the  size  of  digital  mammograms  is  quite  large,  an  ROI  (fixed  at  either  512  x  512  or  1024  x  1024)  within  the  original 
image  was  chosen  to  avoid  computing  over  regions  that  do  not  contain  suspicious  areas.  This  is  also  shown  in  Figure  9, 
where  Figure  9(a)  exhibits  an  original  digitized  mammogram  with  a  512  x  512  ROI  that  contains  a  possible  mass. 


Figure  9(c)  and  Figure  9(d)  display  this  ROI  before  and  after  enhancement  via  non-linear  modification  of  multi-scale 
coefficients,  respectively. 


(a)  (b)  (d) 

Figure  9:  (a)  Original  mammogi^m  with  selected  ROI  containing  a  mass,  (b)  Multi-Scale  Contrest  Enhancement  (MSCE)  GUI, 
(c)  Original  ROI,  and  (d)  Enhanced  ROI. 


3.3.  Display  and  Hardware  Settings 

The  enhancement  protocol  was  executed  on  an  IBM  IntelliStetion  Z  ¥ro  Professional  Workstation  Type  6865.  This  machine 
had  two  Intel  Pentium  n  Xeon  microprocessors  (450  MHz),  512  Mbytes  of  RAM  and  was  equipped  with  36  Gbytes  of  hard 
disk  space.  Windows  NT  4.0  with  service  pack  4  was  the  operating  system. 

To  explore  the  richness  of  information  quanti:red  at  16-bit  per  pixel  (bpp)  grayscale  data  (65536  shades  of  gray),  the  IBM 
IntelliStation  workstation  was  equipped  wifti  two  BARCOMed  5M^1H  Graphics  conteollers.  ITiese  are  hi^-resolution 
display  subsystems  for  the  PCI  bus  with  a  resolution  of 2048x2560  pixels  each,  a  digital-to-analog  converter  (DAC)  capable 
of  1024  shades  of  gray,  and  real  time  window  leveling.  With  the  BARCO  fiamebuflfers,  an  extended  hardware  palette  of 
nearly  16,000  entries  could  be  accessed  through  specialized  ‘‘C”  fiinction  calls  that  were  part  of  a  library  provided  to  us  as 
developers  for  BARCO/Metheus.  Using  these  library  functions,  the  extended  palette  w^  loaded  wifti  a  ramp  of  4096  shades 
of  gray  <x>rresponding  to  12-bit  resolution.  Images  stored  in  16-bit  per  pixel  format,  were  rescaled  to  12  bpp,  if  necessary 
(most  of  the  mammograms  were  digiti^d  at  a  resolution  of  12  bpp),  and  Aen  displayed  at  full  resolution.  Direct  access  to  the 
video  framebuffer  also  sped  up  the  display  process  usefiil  for  uptoing  and  refreshing  the  different  views  on  the  screen. 

Two  high-resolution  MegaScan  monitors  were  attached  to  this  workstetion  providing  dual  heMed  display  on  a  single  logical 
framebuffer  or  virtual  desktop  of  4000x2048  pixels,  respectively  with  Windows  NT  4.0.  To  ensure  iflie  ^curate  depiction  of 
the  same  image  quality  on  both  screens,  a  Bi^CO  P1500  luminance  photometer  was  used.  It  recogni:red  the  1024  shades  of 
gray  displayed  by  a  monitor  and  had  a  range  of  0-450ft-L.  Bofti  monitors  were  calibrated  to  correct  for  non-linearity  of 
display  properties  through  gamma  correction. 

Lighting  conditions  were  controlled  for  the  ROC  study  to  model  reading  room  conditions.  The  ambient  light  intensity  was 
measured  with  the  luminance  photometer  to  be  12.802659  candelea/m^.  It  is  worthwhile  to  note  that  the  optimality  of 
enhancement  pmameters  is  independent  of  the  CRT  display  quality  and  the  image  acquisition  quality.  As  their  computetion  is 
data  driven,  Aey  are  adapted  to  signal  content  and  its  characteristics.  As  our  mdiologiste  gave  us  feedb^k  on  the  quality  of 
die  enhancement,  we  can  adjust  these  initial  default  settings  in  future  studies. 

4.  DESCRIFriON  OF  THE  RECEIVER  OPERATING  CHARACTERISTICS  (ROC)  STUDY 

The  first  receiver  operating  characteristics  (ROC)  study  focused  on  overcomplete  dyadic  wavelets  for  enhancement  of 
mammographic  features  in  digitized  mammograms.  Specifically,  dyadic  spline  wavelet  functions  were  used  together  with  a 
sigmoidal  non-linear  enhancement  ftmction  explicitly  described  in  Section  2.  The  ROC  study  included  three  radiologiste 


specialized  in  mammography.  The  Director  of  the  Breast  Imaging  Center  at  Columbia-Presbyterian  Medical  Center,  Dr. 
Suzanne  Smith,  assisted  in  the  selection  of  cases. 

4.1.  Selection  of  Cases 

To  measure  the  benefits  of  diagnosing  digitized  mammograms  with  enhancement  through  multi-scale  expansions,  we  focused 
on  dense  mammograms,  i.e.  mammograms  of  density  3  and  4  on  the  American  College  of  Radiology  (ACR)  breast  density 
rating,  which  are  the  most  difficult  cases  in  screening.  In  general,  the  enhancement  protocol  aimed  at  improving  the  detection 
and  localization  of  mammographic  features,  such  as  microcalcifications,  masses,  and  spicular  lesions  without  introducing 
“false-positives”. 

To  compare  the  performance  of  radiologists  with  and  without  using  the  enhancement  tool,  two  groups  of  30  cases  each  were 
presented.  Each  group  contained  15  cases  of  cancerous  and  15  cases  of  normal  mammograms.  As  mentioned  above,  a 
national  mammography  database  of  the  University  of  South  Florida  provided  “ground  truth”  (mostly  through  biopsy)  for  the 
selected  cases.  The  selection  was  carried  out  very  carefully  under  the  guidance  of  a  mammographer  (Dr.  Smith),  in  order  to 
find  rather  challenging  cases  of  similar  difficulty  for  each  group.  Images  showing  metal  markers  (“bibis”)  to  indicate 
suspicious  regions  of  breast  tissue  were  avoided  as  well  as  obvious  malignancies.  Due  to  time  constraints  the  number  of  cases 
was  limited  for  this  initial  study. 

4.2.  Paradigm  of  Diagnosis  of  Study 

For  each  case  presented  to  the  radiologist,  the  enhancement  procedure  followed  was  the  following: 

Paradigm  A:  Without  Enhancement: 

The  radiologist  made  a  diagnosis  based  only  on  the  four  original  displays  and  the  viewport.  No  processing  of  ROIs  was 
allowed. 

Paradigm  B:  With  Enhancement: 

The  radiologist  selected  an  ROI  in  one  of  the  views  and  could  apply  multi-scale  enhancement.  Four  levels  of  coefficients 
were  computed.  The  radiologist  then  evaluated  the  quality  of  an  enhanced  ROI  and  adjusted  the  equalizer  sliders  of  a  channel 
to  improve  the  visual  quality  of  suspicious  regions.  Once  he/she  was  satisfied  with  the  visual  result  or  if  he/she  judged  that 
additional  benefit  could  not  be  achieved,  he/she  made  a  diagnostic  decision. 

A  diagnosis  included  specifying  all  lesions  found  and  assigning  a  BI-RAD  scale  to  each  breast  and  the  case.  In  addition,  the 
radiologist  was  asked  to  choose  a  level  of  confidence  (LOG)  for  each  positive  diagnosis,  i.e.  cancer  is  present,  on  an  integer 
scale  from  1  (definitely  negative,  i.e.  total  confidence  that  there  are  no  malignant  lesions)  to  5  (definitely  positive,  i.e.  total 
confidence  that  there  is  a  malignant  lesion).  The  value  for  the  LOG  was  used  in  the  analysis  of  data  to  decide  whether  a  lesion 
was  classified  as  malignant  or  benign  (please  see  discussion  of  LOG  ratings  in  Section  4.4). 

4.3.  ROC  Data 

Table  2  and  Table  3  summarize  the  data  acquired  during  the  study.  Group  1  comprises  the  set  of  cases,  where  the  radiologists 
were  allowed  to  take  advantage  of  the  enhancement  protocol,  whereas  Group  2  contains  those  cases,  where  no  processing 
could  be  applied.  Each  of  the  tables  shows  the  case  numbers,  the  case  designation  and  total  number  (#)  of  lesions  for  each 
case  according  to  the  mammography  database  (DB),  and  for  each  of  the  three  mammographers  the  BI  RAD  rating  and  level 
of  confidence  (LOG)  values.  The  BI_RAD  rating  could  be  chosen  from  the  standard  categories  0-5  with  0  meaning  that 
additional  information  for  a  more  confident  diagnosis  was  needed.  In  such  cases,  the  radiologists  were  asked  to  also  select  a 
BI  RAD  rating  different  from  0,  if  they  were  asked  to  make  a  diagnosis  without  any  additional  information.  This  number  is 
shown  in  parentheses  for  such  cases. 

In  each  table  both  groups  are  sorted  into  actually-negative  cases  (normals  with  “0  lesions)  and  actually-positive  cases 
(cancers  with,  at  least  “1  malignant  lesion),  since  this  is  required  for  subsequent  analysis  of  the  data. 


Table  2:  ROC  data  for  three  mammographeis  for  Group  1,  Le.  with  Enhancement  enabled. 
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Table  3:  ROC  data  for  three  mammogi^phei^  for  Group  2,  Le.  without  enhancement. 


4.4.  ROC  Analysis:  General  Principles 

The  most  widely  used  method  to  objectively  evaluate  the  performance  of  a  diagnostic  system  or  the  difference  in 
performance  between  two  diagnostic  systems  is  ROC  analysis.  It  compares  radiologists’  image-based  diagnoses  with  known 
states  of  disease  and  health.  In  ROC  analysis,  performance  of  a  diagnostic  system  is  described  by  the  indices  of  “sensitivity” 
and  “specificity”,  where  “sensitivity”  can  be  expressed  as  the  true-positive  fraction  (TPF)  and  “specificity”  by  the  true¬ 
negative  fraction  (TNF)  of  a  diagnosis  [16].  In  a  complimentary  way,  the  false-negative  fraction  (FNF)  and  the  false-positive 
fraction  (FPF)  can  be  defined  as  FNF  =  1-TPF  and  FPF  =  1-TPF,  respectively,  with  a  similar  interpretation.  Due  to  this 
dependence,  it  is  only  necessary  to  measure  one  pair  of  indices,  and  frequently  TPF  and  FPF  are  used  (as  in  our  study). 

The  underlying  model  for  ROC  analysis  is  the  use  of  probability  density  distributions  of  a  radiologist’s  confidence  in  a 
positive  diagnosis  for  a  particular  diagnostic  task  for  true  positive  and  true  negative  patients  [16].  It  is  currently  accepted  that 
based  on  a  confidence  threshold,  i.e.  a  particular  level  of  confidence  (LOC)  in  a  positive  diagnosis,  a  diagnosis  is  considered 
to  be  positive,  if  it  exceeds  this  threshold,  and  a  diagnosis  is  considered  to  be  negative,  if  it  falls  below  the  threshold.  TPF  and 
FPF  are  then  calculated  from  the  probability  density  distributions  as  areas  under  the  curves  delimited  by  the  confidence 
threshold  (see  Figure  10).  If  the  confidence  threshold  is  varied  continuously,  an  ROC  curve  can  be  generated  from  the  pair 
values  for  TPF  and  FPF.  ROC  curves  that  indicate  better  decision  performance  are  positioned  higher  in  the  unit  square 
spanned  by  FPF  and  TPF  (higher  TPF  values  for  the  same  FPF  values).  The  area  under  the  ROC  curve,  provides  a  useful 
summary  index  for  the  inherent  discrimination  performance  of  a  diagnostic  system.  Thus,  A2  is  the  average  value  of 
sensitivity  of  a  corresponding  ROC  curve,  if  the  specificity  of  the  system  is  selected  randomly  between  0.0  and  1.0. 
Conversely,  it  can  be  considered  as  the  average  value  specificity  of  a  corresponding  ROC  curve,  if  the  sensitivity  of  the 
system  is  selected  randomly  between  0.0  and  1.0  [16]. 
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Figure  10:  Schematic  example  of  the  model  that  underlies  ROC  analysis.  The  bell-shaped  curves  represent  probability  density 
distributions  of  a  radiologist's  confidence  in  a  positive  diagnosis.  A  confidence  threshold,  represented  by  a  vertical  line, 
separates  “positive”  decisions  from  “negative”  decisions  (This  figure  was  reprinted  from  [16]). 

In  practice,  data  for  an  ROC  analysis  is  obtained  by  providing  a  set  of  rating  categories  to  the  radiologist.  For  a  rating  scale 
we  chose  discrete  values  from  1  to  5  for  the  level  of  confidence  (LOC)  in  a  positive  diagnosis.  The  meaning  of  these  values 
was  as  follows:  (1)  definitely  or  almost  definitely  negative,  (2)  probably  negative,  (3)  possibly  positive,  (4)  probably  positive, 
and  (5)  definitely  or  almost  definitely  positive.  With  this  choice  the  value  for  the  LOC  is  similar  to  the  standard  BI  RAD 
rating  scale  used  in  screening. 

To  generate  an  ROC  curve  from  discrete  data  requires  assumptions  about  the  functional  form  of  the  curve.  The  “binormal” 
model  has  been  widely  used  in  medical  imaging.  This  model  includes  two  adjustable  parameters,  and  it  is  assumed  that  each 
conventional  ROC  curve  has  the  same  functional  form  as  that  implied  by  two  “normal”  (i.e.,  Gaussian)  decision  variable 
distributions  with  generally  different  means  and  standard  deviations  [37],  [38]. 

The  two  adjustable  parameters  of  the  binormal  ROC  curve  can  be  taken  to  be  the  y-intercept  and  the  slope  of  the  straight  line 
that  represents  the  ROC  curve,  when  it  is  plotted  on  normal-deviate  axes.  These  two  parameters,  denoted  as  “a”  and  “Z>”,  can 
be  interpreted  as  an  effective  pair  of  underlying  Gaussian  distributions  as  the  distance  between  the  means  of  the  two 
distributions  and  the  standard  deviation  of  the  actually  negative  distribution,  respectively  with  both  expressed  in  units  of  the 
standard  deviation  of  the  actually  positive  distribution  [16].  With  the  binormal  model,  a  maximum-likelihood  parameter 
estimation  scheme  is  then  used  to  generate  an  ROC  curve  that  best  represents  the  data. 

If  two  different  diagnostic  systems  are  to  be  evaluated,  the  statistical  difference  of  an  apparent  difference  between  measured 
ROC  curves  is  of  interest.  Testing  differences  between  ROC  curves  is  well  described  in  the  literature  [39],  [40]. 


4.S.  Results  from  ROC  Anal^is 

In  our  study,  ROC  analysis  was  possible,  since  the  ‘‘ground  truth”  for  each  case  was  provided  by  the  mammogmphy  datobase. 
In  general,  any  enhancement  protocol  should  incre^e  sensitivity,  Le.  fraction  of  true-positives  (TPF),  without  decreasing 
specificity,  i.e.  essentially  without  incroasing  the  fi^ction  of  false-positives  (FPF)  [41].  An  initial  analysis  of  the  data  counted 
the  number  of  false-positives  and  true-positives  in  each  group  of  cases.  Before  a  lesion  was  considered  being  diagnosed  as 
malignant  or  benign,  the  LOC  value  was  thresholded  [16].  The  fiireshold  value  influence  the  shape  of  the  ROC  curve  and  its 
interpretation.  For  example,  if  the  fiireshold  for  the  level  of  confidence  was  chosen  to  be  3,  meaning  that  lesions  with  a  LOC 
greater  or  equal  3  were  cronsiderad  as  malignant  then  the  average  TPF  was  found  to  be  0.667  wifii  enhancement,  and  TPF  = 
0.569  without  enhancement,  Tliis  observed  increase  in  sensitivity  is  encouraging,  though  it  was  accompanied  by  a  slight 
increase  in  the  fraction  of  false-positives  (0.222  compared  to  0.17S).  The  latter  is  not  too  surprising,  since  the  <^plied 
enhancement  protocol  only  used  dyadic  spline  wavelets  with  the  non-linear  sigmoidal  enhancement  fiinction,  which  is 
certainly  not  optimal  for  il  types  of  lesions.  We  believe  that  dyadic  spline  wavelet  expansions  are  best  used  to  enhance 
microcalcifications.  If  the  analysis  of  the  data  only  focused  on  microcalcifications,  fiien  we  observed  TPF  =  0.417  with 
enhmicement  compared  to  TPF  =  0.222  without  enhancement.  No  increase  or  decre^e  in  FPF  was  noticed!  The  last  finding 
supports  the  promise  for  future  research  to  design  specific  enhancement  protocols  for  each  mammographic  feature.  Table  4 
summaries  initial  results  of  toe  ROC  study  using  the  single  basis  function  described  earlier  in  Section  2.3. 


Table  4:  Results  of  prelimlnaiy  ROC  study.  TPF  refera  to  the  fraction  of  true-positives  and  FPF  to  the  fi'action  of  false-positives. 


A  more  thorough  analysis  of  the  data  was  undertoken  by  using  the  ROCKIT  software  developed  by  a  research  group  led  by 
Charles  Metz  at  the  Univeraity  of  Chicago  [42],  [43].  This  software  package  was  written  to  analyre  data  from  ROC  studies 
and  to  generate  corresponding  ROC  curves.  More  specifically,  the  purpose  of  ROCKIT  is  to  calculate  maximum-likelihood 
estimates  of  the  parameter  of  a  conventional  “binomial”  model  for  the  input  data,  to  calculate  maximum-likelihood 
estimates  of  the  parameters  of  a  “bivariate  binormal”  model  for  data  from  two  potentially  correlated  diagnostic  teste  and, 
thus,  to  estimate  the  binormal  ROC  curves  implied  by  those  data  and  their  itorrelation;  and  to  calculate  the  stotistical 
significance  of  the  difference  between  two  ROC  curve  estimates  using  any  one  of  three  distinct  statistical  teste: 

1.  The  Bivamte  Test:  A  bivariate  Chi-square  test  of  the  simultaneous  differences  between  the  “a”  parameter  and 
between  file  “5”  parameters  of  the  two  ROC  curves.  {Null  hypothesis:  the  data  sets  arose  from  the  same  binormal 
ROC  curve.) 

2.  The  Area  Test:  A  univariate  z-score  test  of  the  difference  between  the  areas  under  the  two  ROC  curves.  {Nutt 
liypotliesis:  the  data  sets  arose  tom  binormal  ROC  curves  with  equal  m*eas  beneath  them.) 

3.  The  TFP  Test:  A  univariate  z-scorc  test  of  the  difference  between  the  true-positive  fractions  (TPFs)  on  the  two  ROC 
curves  at  a  selected  false-positive  fi^ion  (FPF).  (Null  hypothesis:  the  data  sets  arose  from  binormal  ROC  curves 
having  the  same  TPF  at  tiie  selected  FPF.) 

Three  types  of  input  data  are  allowed  for  stotistical  testing  of  file  differences  between  ROC  curves: 

1.  Unpaired  (uncorrelated)  test  resulte.  The  two  “cxmditions”  are  applied  to  independent  case  samples  — ^  for  example, 
tom  two  different  diagnostic  tests  performed  on  fiie  different  patiente,  tom  two  different  i^iologists  who  make 
probability  judgments  concerning  the  presence  of  a  specified  disease  in  different  images,  etc.; 

2.  Fully  paired  (correlated)  test  results,  in  which  data  from  botfi  of  two  conditions  are  available  for  each  case  in  a  single 
case  sample.  The  two  “<x>nditions”  in  each  test-result  pair  could  correspond,  for  e^rample,  to  two  different  di^ostic 
tests  performed  on  the  same  patient,  to  two  different  radiologists  who  make  probability  judgmente  concerning  the 
presence  of  a  specified  dise^e  in  file  same  image,  etc.;  and 

3.  Partially-paired  test  results  —  for  example,  two  different  diagnostic  teste  performed  on  flie  same  patient  sample  and 
on  some  additional  patients  who  r^eived  only  one  of  the  diagnostic  tests. 


Rockit  assumes  that  the  population  ROC  curve  for  each  condition  plots  as  a  straight  line  on  “normal-deviate”  axes,  or 
equivalently,  that  the  input  data  follow  normal  distributions  after  some  unknown  monotonic  transformation  [16].  ROC  curves 
measured  in  a  broad  variety  of  fields  demonstrate  this  “binormal”  form  [44],  [45],  and  [46].  The  assumption  may  be  satisfied 
even  when  the  raw  data  have  multimodal  and/or  skewed  distributions  [43],  [42]. 

Using  the  ROCKIT  software  the  analysis  was  first  applied  independently  to  the  datasets  for  Group  1  and  Group  2  for  each  of 
the  three  radiologists.  Unfortunately,  this  approach  did  not  allow  us  to  compare  the  diagnostic  performance  for  the  two 
diagnostic  systems  (softcopy  display  with  and  without  enhancement).  The  reason  for  that  was  that  the  analysis  for,  at  least 
one  group  of  cases  could  not  be  completed,  since  the  data  was  found  to  be  degenerate  [41].  In  this  case,  the  result  of  the  ROC 
analysis  would  be  a  straight  line  with  a  constant  value  for  TPF,  and,  therefore  the  software  aborts  processing  to  avoid 
meaningless  output.  According  to  the  authors  of  the  software,  a  degenerate  data  distribution  can  be  found,  if  the  number  of 
samples  is  too  small  or  in  datasets  with  many  tied  values  [43]. 

Since  the  number  of  cases  could  not  be  increased  after  conducting  the  study,  and  in  order  to  obtain  more  complete  results,  we 
decided  to  apply  the  analysis  to  the  union  of  data  from  all  three  radiologists.  This  was  justified  by  the  fact  that  all  three 
radiologists  came  fi-om  the  same  population  with  a  similar  level  of  experience.  Thus,  their  performance  should  be  similar 
under  the  same  conditions,  and  the  data  could  be  treated  as  independent  samples  (unpaired  data).  If  the  data  did  not  have  to 
be  pooled,  it  would  have  been  unpaired,  since  the  two  different  conditions  were  applied  to  different  sample  cases. 
Nevertheless,  we  are  well  aware  that  the  statistical  significance  of  the  results  must  be  interpreted  carefully.  For  future  ROC 
studies  we  plan  to  increase  the  number  of  cases,  in  order  to  avoid  such  a  problem.  To  check  on  our  assumption  of 
independent  samples  (unpaired  data)  and  for  completion  we  also  repeated  the  analysis  with  the  input  as  paired  data.  These 
results  are  included  in  this  chapter  as  well. 

For  the  analysis  Group  1  (with  enhancement)  was  set  as  Condition  1  and  Group  2  (without  enhancement)  was  considered  as 
Condition  2.  The  resulting  ROC  curves  for  data  analyzed  as  unpaired  are  shown  in  IPigure  11.  Their  corresponding  values  for 
FPF  and  TPF  are  given  in  Table  5.  Finally,  the  most  important  results  of  ROC  analysis,  the  binormal  parameters  a,  b,  and  the 
area  under  the  ROC  curve  with  their  corresponding  standard  errors,  95%  confidence  intervals,  and  correlation  of  a  and  b 
are  summarized  for  unpaired  data  in  Table  6.  Note  that  the  95%  confidence  intervals  are  symmetric  for  the  binormal 
parameters  a  and  b,  but  asymmetric  for  the  area  index  A2.  The  corresponding  results  firom  the  analysis  as  paired  data  follow 
directly  afterwards.  ROC  curves  are  shown  in  Figure  12,  FPF  and  TPF  values  in  Table  7,  and  parameters  a,  b,  and  A^  together 
with  their  corresponding  standard  errors,  95%  confidence  intervals,  and  correlation  of  a  and  b  in  Table  8. 

ROC  Curves  for  Data  with  and  without  Enhancement 


False  Positive  Fraction  (FPF) 


— ^With  Enhancement 
-^Without  Enhancement 


Figure  11:  ROC  curves  for  data  with  Condition  1  (with  enhancement)  and  Condition  2  (without  enhancement)  analyzed  as 
unpaired  data  (independent  analysis). 


Table  5t  Values  for  fiilse-positive  fractions  (FPE)  and  true-positive  fractions  (TPF)  for  Condition  1  (with  enhancement,  TPF  1)  and 
Condition  2  (without  enhancement  TPF  2)  ana^ed  as  unpaired  date  (independent  anaty^sis). 
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Figure  12:  ROC  curves  for  data  with  Condition  1  (with  enhancement)  and  Condition  2  (without  enhancement)  analyzed  as  paired 
data  (correlated  analysis). 
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Table  7:  Values  for  false-positive  fractions  (FPF)  and  true-positive  fractions  (TPF)  for  Condition  1  (with  enhancement,  TPF  1)  and 
Condition  2  (without  enhancement,  TPF  2)  analyzed  as  paired  data  (correlated  analysis). 
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Table  Binonnal  parameters  b,  area  under  ROC  curve  with  their  cori^sponding  standard  erroi^,  95%  confidence  intervals, 
and  con^lation(a,  b)  for  Condition  1  (with  enhancement)  and  Condition  2  (without  enhancement)  analyzed  as  paired  data 
(correlated  analysis), 

4.6.  Discussion 

As  seen  fiom  the  analysis  for  unpaired  data,  the  value  for  the  area  under  the  ROC  curve  was  by  8.7%  larger  for  Q>ndition 
1  (wilh  enhancement)  than  it  for  Condition  2  (without  enhancement).  In  all  c^es  the  standard  eitor  for^^  was  between 
0.03  md  0.05,  which  w^  rather  small.  Though  the  95%  confidence  intervals  for  overlapped,  there  was  a  clear  tendency 
that  diagnostic  performance  improved  with  enhancement  in  comparison  wifii  diagnosis  without  enhancement  All  ROC 
curves  lay  high  in  the  unit  square  of  FPF  and  TPF,  which  corresponded  to  accurate  diagnostic  performances  in  general,  but 
the  curve  for  Cbndition  1  was  positioned  slightly  higher  (see  Figure  11). 

Similar  resulte  were  generally  obtoined  for  the  analysis  as  paired  data.  The  increase  in  for  Condition  1  with  respect  to 
Condition  2  was  8.5  %,  but  fiiere  was  an  overlap  of  the  95%  confidence  intervals  for  A^  as  well.  The  ROC  curve  for 
Condition  1  was  also  positioned  slightly  higher  than  the  one  for  Condition  2  (see  Figure  12).  Values  for  a,  h,  and  A^  were 
very  similar  for  both  types  of  analysis.  Hence,  the  same  tendency  of  improved  di^ostic  performance  with  enhancement 
compared  to  dia^osis  without  can  be  inferred. 

The  observed  increase  of  the  summmy  index  A^  within  statistical  errors  and  the  higher  position  of  fiie  ROC  curve  for 
diagnosis  with  enhancement  encourage  us  to  ftirtfier  puraue  the  application  of  enhancement  protocols  for  mammogiaphic 
screening.  We  are  aware  of  the  fact  that  fiiere  always  are  inherent  sources  of  variability  in  the  index  A^,  such  m  a  ‘^case- 
sample’’  <x>mponent  due  to  random  variations  in  the  difficulty  of  the  cases  included  in  an  ROC  experiment,  a  “between- 
reader”  component  due  to  random  variations  in  the  skills  of  the  observers  participating  in  the  experiment,  and  a  ‘%ithin- 
reader”  component  associated  with  each  reader’s  inability  to  reproduce  her^is  diagnosis  of  eveiy  case  on  repeated  readings 
[16].  hi  addition,  we  were  not  able  to  analyze  the  data  for  each  radiologist  separately  due  to  date  degeneracy  as  mentioned 
above.  The  latter  has  diminished  the  statistical  significance  of  our  results  obtained  from  die  analysis  of  all  data  combined, 
since  not  all  samples  were  completely  independent. 

Hence,  for  firture  ROC  studies  we  plan  to  increase  the  number  of  c^es  to  avoid  degenerate  dateets  for  the  analysis  and  to 
increase  Ae  statistical  power  of  the  experiment. 

Aside  from  statistical  considerations  and  Ae  cautious  inteipretation  of  the  results  of  Ais  study  we  know  Aat  our  prototype 
test  bed  software  tool  can  be  fiirther  optimized.  To  improve  multi-scale  contot  enhancement  Ae  idea  is  to  develop  feature 
specific  enhancement  protocols  with  different  bases  and  ^sociated  non-linear  functions  for  each  distinct  mammographic 
feature,  such  as  microcalcifications,  masses,  and  spicular  lesions.  The  enhancement  protocol  used  for  Ais  experiment,  dyadic 
spline  wavelets  wiA  non-lmear  sigmoidal  fimction,  w^  suggested  to  work  best  for  microcalcifications  according  to  our 
previous  work  with  multi-scale  expansions  [2],  [25].  The  results  of  Ais  first  ROC  experiment  seem  to  confirm  our 
expectations. 


5,  CONCLUSIONS  AND  FUTURE  WORK 

We  have  reported  on  the  successful  completion  of  the  first  receiver  operating  characteristics  (ROC)  study  to  evaluate  the 
benefits  of  contrast  enhancement  via  overcomplete  multi-scale  expansions  of  mammograms.  The  study  was  carried  out  in 
collaboration  with  radiologists  at  the  Breast  Imaging  Center  in  Columbia-Presbyterian  Medical  Center  and  the  Biomedical 
Imaging  Laboratory  of  Columbia  University. 

In  continuation  of  our  previous  work  in  digital  mammography,  an  enhancement  protocol  using  a  dyadic  spline  wavelet  as  the 
basis  for  multi-scale  expansion  and  an  associated  non-linear  sigmoidal  enhancement  function  was  designed.  Suspicious  areas 
(ROIs)  of  digitized  mammograms  were  decomposed  onto  a  multi-scale  basis  to  obtain  coefficients  at  distinct  subbands. 
Coefficients  were  modified  by  applying  a  non-linear  sigmoidal  function.  Two  parameters  could  be  adjusted  to  change  the 
nature  of  enhancement.  Image  reconstruction  from  modified  coefficients  occurred  in  nearly  real  time  through  an  interactive 
interface  running  on  a  high-resolution  digital  mammography  workstation.  To  visualize  raw  data  of  digitized  mammograms  at 
the  highest  possible  contrast  and  spatial  resolutions,  16-Bit  BARCO/Metheus  framebuffers  together  with  a  dual  headed  high- 
resolution  MegaScan  grayscale  monitor  were  utilized  in  hardware.  We  incorporated  specialized  software  function  calls  to 
directly  access  the  video  framebuffer  for  fast/smooth  image  display  and  update. 

To  quantify  the  performance  of  our  multi-scale  based  processing  technique  in  terms  of  overall  sensitivity  and  specificity,  an 
ROC  study  was  designed  and  conducted  with  three  radiologists  from  Columbia-Presbyterian  Medical  Center  specialized  in 
mammography.  Conventional  ROC  curves  were  generated  and  significant  statistical  parameters  determined.  The  area  under 
the  ROC  curve  was  used  as  a  summary  index  to  quantify  overall  specificity  and  sensitivity  of  the  two  diagnostic  systems 
[16].  Unfortunately,  it  was  not  possible  to  analyze  datasets  for  each  of  three  mammographers  separately  due  to  data 
degeneracy.  Nevertheless,  analyzing  all  the  data  together  yielded  a  slight  increase  (8.7%)  in  the  area-^^  for  diagnosis  with 
enhancement  compared  to  diagnosis  without.  Despite  the  limited  statistical  significance  of  this  result,  it  encourages  us  to 
further  investigate  the  application  of  multi-scale  methods  for  contrast  enhancement  of  mammograms.  More  extensive  ROC 
studies  with  a  larger  number  of  cases  are  planned  to  further  evaluate  the  benefits  of  such  processing  techniques. 

Ancillary  to  statistical  results,  we  received  very  positive  feedback  from  the  participating  radiologists,  who  expressed  great 
interest  in  using  the  interactive  display  tool  and  acknowledged  a  marked  improvement  in  image  quality,  when  enhancement 
was  applied. 

The  current  enhancement  protocol  works  best  for  the  detection/enhancement  of  microcalcifications.  Future  directions  of  work 
include  the  expansion  of  the  choice  of  enhancement  protocols  to  a  menu  of  feature  specific  enhancement  algorithms  tailored 
for  each  mammographic  feature,  such  as  microcalcifications,  masses,  and  spicular  lesions,  e.g.  the  application  of  brushlet 
functions  [47],  [48]  to  mammograms  with  spicular  lesions.  In  addition,  the  investigation  of  a  range  of  optimal  enhancement 
parameters  and  the  optimization  of  our  interface  software  tool  comprise  further  projects.  Our  “dream”  is  to  present  a  clinical 
interface,  where  specific  enhancement  protocols  can  be  selected  by  a  physician  by  only  “pushing  a  button  on  the  screen”.  We 
envision  that  through  such  a  clinical  interface  the  diagnostic  performance  of  radiologists  in  screening  digital  mammograms 
could  be  substantially  improved,  both  in  terms  of  cost  and  quality. 
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Image  enhancement  in  mammography  is  typically  concerned  with  either  general  visibility 
of  all  features  or  conspicuity  of  a  specific  sign  of  malignancy.  We  describe  a  synthesis  of 
the  two  approaches  through  fusion  of  locally  enhanced  microcalcifications,  circumscribed 
masses,  and  stellate  lesions.  Both  local  processing  and  image  fusion  are  performed  within 
a  single  wavelet  transform  framework  which  contributes  to  the  computational  efficiency 
of  the  method.  The  algorithm  not  only  allows  for  efficient  combination  of  specific  fea¬ 
tures  of  importance,  but  also  provides  a  flexible  framework  for  incorporation  of  distinct 
enhancement  methods  and  their  independent  optimization. 

Mammography,  contrast  enhancement,  image  fusion,  wavelet  transform. 


1.  INTRODUCTION 

In  general,  mammographic  image  enhancement  methods  target  either  visualization  of 
all  features  in  an  image  [1,  2,  3,  4]  or  visibility  of  specifle  features  of  importance  such  as 
microcalcifications  [5]. 

Methods  from  the  first  category  are  not  optimized  for  a  specific  type  of  cancer,  and 
are  often  developed  for  a  framework  more  general  than  mammography  alone.  The  second 
category  approaches  can  be  quite  successful  in  their  area  of  specialization;  however,  in 
order  to  process  mammograms  for  presence  of  distinct  features,  independent  application 
of  different  algorithms  could  result  in  both  larger  number  of  images  to  be  interpreted  by 
a  radiologist  and  increased  computational  complexity. 

Here,  we  present  an  approach  which  overcomes  these  shortcomings  and  problematic 
limitations  via  synthesis  of  the  two  paradigms  by  means  of  image  fusion.  The  algorithm 
consists  of  two  major  steps:  (1)  wavelet  coefficients  are  modified  distinctly  for  each  type 
of  malignancy;  (2)  the  obtained  multiple  sets  of  wavelet  coefficients  are  fused  into  a  single 
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set  from  which  a  reconstruction  is  computed.  The  scheme  allows  efficient  deployment  of 
an  enhancement  strategy  appropriate  for  clinical  screening  protocols:  an  enhancement 
algorithm  is  first  developed  for  each  specific  type  of  feature  independently,  and  the  results 
are  then  combined  using  an  appropriate  fusion  strategy. 

2.  WAVELET  TRANSFORM 


Wavelet  based  methods  are  particularly  well  suited  for  processing  of  mammograms  since 
mammographic  features  greatly  vary  in  shape  and  size.  Commonly  used  orthogonal  and 
biorthogonal  wavelet  transforms,  however,  may  not  be  the  best  tool  for  mammographic 
image  enhancement  because  their  lack  of  translation  invariance  can  lead  to  artifacts  possi¬ 
bly  affecting  a  radiologist’s  interpretation.  Translation-invariant  but  overcomplete  wavelet 
representations  avoid  artifacts  and  have  been  successfully  used  for  processing  of  mammo¬ 
grams  [1,  2,  5]. 

Rotation  invariance  is  another  desirable  property  of  wavelet  decompositions.  The 
concept  of  steerability  [6]  has  been  utilized  for  construction  of  wavelet  transforms  enabling 
rotation-invariant  processing  of  mammograms  [3].  Our  scheme  is  built  around  a  multiscale 
spline  derivative-based  transform  which,  in  addition  to  being  translation-invariant  and 
approximately  steerable,  is  also  suitable  for  non-linear  methods  of  enhancement. 

We  use  x-y  separable  wavelets 


where  0p{x)  denotes  a  central  B-spline  of  order  p,  and  limit  ourselves  to  first  and  second 
derivatives  d  G  {1,2}.  Figure  1  shows  wavelets  with  p=3. 

A  rotation  of  wavelet  by  angle  d  can  be  expressed  as 


i^ix^y) 


d-ii  ^Pp+d{x)  (f^p+djy) 


n“  n 


dx^~^ 


where  ft  =  {(X3S0,sm0)  =  {nx,ny).  The  terms  ^  represent  basis  functions 

needed  to  approximately  steer  wavelet  ilf{x,y).  A  dyadic  wavelet  transform  using  these 
basis  functions  can  be  implemented  as  a  filter  bank  consisting  of  one-dimensional  filters 
only  [7]. 


3.  ENHANCEMENT  OP  MAMMOGRAPHIC  FEATURES 
3.1.  Microcalcifications 

Microcalcifications  appear  on  mammograms  in  approximately  half  of  breast  cancer 
cases.  The  assessment  of  shape,  number,  and  distribution  of  microcalcifications  is  impor¬ 
tant  for  a  rMiologist  to  reach  diagnosis.  Microcalcifications  are  smaller  than  1  mm  in  size 
and  can  be  difficult  to  locate  when  they  are  superimposed  on  dense  breast  tissue. 

Several  techniques  have  been  developed  to  improve  the  visibility  of  microcalcifications 
[5,  8,  9].  The  approach  devised  by  Strickland  and  Hahn  [5]  is  particularly  well  suited 


Figure  1.  Spline  derivatives  in  the  a;-axis  direction,  (a)  Wavelet  equal  to  the  first 
derivative  of  a  quartic  spline,  (b)  Wavelet  equal  to  the  second  derivative  of  a  quintic 
spline. 


for  our  framework:  they  used  an  undecimated  wavelet  transform  to  approximate  second 
derivatives  of  a  Gaussian  probability  density  function  for  a  multiscale  matched  filtering 
for  presence  of  microcalcifications. 

Strickland  and  Hahn  based  their  method  on  the  observation  that  the  average  mi¬ 
crocalcification  can  be  modeled  by  a  circularly  symmetric  Gaussian  function.  We  take 
advantage  of  this  fact  to  model  microcalcifications  by  central  B-splines.  Using  B-spline 
approximations  of  a  Gaussian  function,  the  assumption  that  a  Gaussian  object  is  visible 
approximately  over  ±(7  pixels  [5],  and  the  fact  that  mammograms  in  the  University  of 
Florida  database  were  digitized  at  116/xm  resolution,  four  levels  of  the  transform  described 
in  Section  2  with,  for  example,  p=3  are  needed  to  encompass  different  sizes  of  microcal¬ 
cifications.  The  wavelet  decomposition  including  voices  at  scales  3  and  6  (corresponding 
to  Strickland  and  Hahn’s  octaves  “2.5”  and  “3.5”)  was  obtained  from  relations  between 
central  B-splines  at  integer  scales  [10]. 

The  wavelet  decomposition  enables  approximations  both  to  the  second  derivatives  of 
Gaussian  along  x  and  y  directions  and  to  Laplacian  of  Gaussian  across  distinct  scales 
employed  by  Strickland  and  Hahn.  We  proceed  in  a  similar  fashion:  the  two  outputs 
per  scale  are  thresholded  independently,  all  binary  results  are  then  combined,  a  circular 
region  centered  at  detected  pixel  locations  is  next  multiplied  by  a  gain,  and,  finally,  the 
modified  transform  coefficients  are  used  for  image  fusion. 

3.2.  Circumscribed  Masses 

Almost  half  of  missed  cancers  appear  on  mammograms  as  masses.  Perception  is  a  problem 
particularly  for  patients  with  dense  fibroglandular  patterns.  The  detection  of  masses  can 
be  especially  difficult  because  of  their  small  size  and  subtle  contrast  compared  with  normal 
breast  structures. 


Fan  and  Laine  [2]  developed  a  discrete  dyadic  wavelet  transform  based  algorithm 
suitable  for  enhancement  of  masses.  They  constructed  an  approximation  to  Laplacian  of 
Gaussian  across  dyadic  scales  for  an  isotropic  input  to  a  piecewise  linear  enhancement 
fiinction. 

An  approximation  to  a  Laplacian  of  Gaussian  across  dyadic  scales  te  easy  to  obtain 
using  multiscale  spline  derivatives  from  Section  2:  basis  functions  ^ 
d  =  2  and  i  €  {1,2}  approximate  the  second  derivative  of  a  Gaussian  function  along 
directions  of  x  and  y  axis.  The  appropriate  transform  coefficient  at  each  dyadic  scales  are 
then  added  and  their  sum  input  to  the  piecewise  linear  function  [2] 


X  -  {K-1)T  if  X  <  -T 
Cix)  =  {  Kx  if  |a;l  <  T 

a;  +  iK-l)T  itx  >  T 


used  at  each  level  m+1  of  the  transform  separately.  Due  to  the  expected  size  of  masses, 
levels  greater  than  4  are  enhanced  more  aggressively. 

The  multiplicative  factor  obtained  as  the  ratio  between  the  output  and  input  of  the 
enhancement  function  is  next  applied  to  the  original  wavelet  coefficients  before  fusion  and 
the  associated  inverse  wavelet  transform  are  carried  out. 


3.3.  Stellate  Lesions 

It  is  important  for  radiologists  to  identify  stellate  lesions  since  their  presence  is  a  serious 
indicator  of  malignancy.  Stellate  lesions  vary  in  size  and  subtlety  and,  in  addition,  do  not 
have  a  clear  boundary,  making  them  difficult  to  detect. 

In  the  development  of  our  algorithm,  we  utilized  an  observation  made  by  Kegelmeyer 
et  al.  about  the  distortion  of  edge  orientation  distribution  induced  by  a  stellate  lesion  [11]. 
Normal  mammograms  show  a  roughly  radial  pattern  with  structure  radiating  from  the 
nipple  to  the  chest  wall.  A  stellate  lesion  not  only  changes  this  pattern,  but  also  creates 
another  center  from  which  rays  radiate. 

The  wavelet  transform  from  Section  2  allows  directional  analysis  using  approximations 
to  both  first  and  second  steerable  derivatives  of  a  Gaussian.  A  multiscale  derivative-pair 
quadratic  feature  detector  was  computed  by  finding  the  maximum  of  the  local  oriented 
energy  with  respect  to  angle  0, 

Elm{x,y)  =  sJ{Wl%ns{x,y)Y  +  {W2^ms{x,y)f,  (2) 

where  Wllms{x,y)  and  W22ms(x,  ?/)  denote  wavelet  decompositions  using  first  (Equation 
(1)  with  d=l)  and  second  (Equation  (1)  with  d  =  2)  derivati^  wavelet,  respectively, 
steered  to  angle  9.  The  angle  that  maximizes  the  local  oriented  energy  (2)  represents 
orientation  at  pixel  location  (a:,  y). 

Similar  to  the  method  from  Section  3.1,  processing  is  carried  out  within  windows  of 
scale  dependent  size:  1-norm  of  differences  between  the  local  and  average  orientations 
was  computed  in  the  window  and  used  as  a  measure  of  orientation  nonuniformity.  Soft 
thresholding  as  a  fiinction  of  the  orientation  nonuniformity  measure  was  next  applied  to 
the  transform  coefficients  at  each  dyadic  scale  independently.  The  altered  coefficients  are 
then  included  for  fusion  and  reconstruction. 


Figure  2.  (a)  The  cranio-caudal  view  of  the  left  breast,  (b)  Enhanced  image  improves 
visualization  of  the  borders  of  the  mass. 


4.  FUSION  OF  ENHANCED  FEATURES 

After  coefficients  are  processed  for  enhancement  of  distinct  mammographic  features,  the 
corresponding  coefficients  are  combined  according  to  a  fusion  rule  into  a  new  set  of  trans¬ 
form  coefficients  from  which  the  fused  result  is  reconstructed.  As  a  fusion  rule,  the  max¬ 
imum  oriented  energy  criterion  was  chosen:  at  each  position  and  scale  of  the  transform, 
the  coefficient  with  greatest  local  energy  was  selected  [12]. 

It  is  also  possible  to  put  distinct  weights  on  selected  features,  and  exclude  other 
features  from  the  final  result. 

Figure  2  shows  the  original  mammogram  and  the  processed  image  with  improved  con¬ 
trast  between  the  fat  and  glandular  tissue. 


5.  CONCLUSION 

The  presented  method  incorporates  a  variety  of  properties  of  mammographic  image  en¬ 
hancement  techniques  tailored  to  specific  signs  of  malignancy  into  a  unified  computa¬ 
tional  framework.  A  multiscale  spline  derivative-based  transform  proved  flexible  enough 
for  implicit  enhancement  of  individual  types  of  mammographic  features  and  thus  enabled 
processing  within  a  single  wavelet  transform  decomposition.  In  addition  to  its  efficiency, 
the  algorithm  is  well  suited  for  further  refinements;  optimizations  can  be  performed  for 
each  type  of  malignancy  alone,  and  separately  for  the  fusion  strategy. 

Our  preliminary  experiments  imply  that  an  enhancement  via  fusion  approach  can  pro¬ 
vide  more  obvious  clues  for  radiologists.  Further  clinical  tests  are  planned  to  verify  that 
the  versatility  of  this  paradigm  can  provide  a  better  viewing  environment  for  a  more  re¬ 
liable  interpretation  in  screening  mammography. 
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